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Abstract. This paper concerns with mesh restrictions that are needed to satisfy several impor¬ 
tant mathematical properties - maximum principles, comparison principles, and the non-negative 
constraint - for a general linear second-order elliptic partial differential equation. We critically 
review some recent developments in the field of discrete maximum principles, derive new results, 
and discuss some possible future research directions in this area. In particular, we derive restric¬ 
tions for a three-node triangular (T3) element and a four-node quadrilateral (Q4) element to satisfy 
comparison principles, maximum principles, and the non-negative constraint under the standard 
single-field Galerkin formulation. Analysis is restricted to uniformly elliptic linear differential op¬ 
erators in divergence form with Dirichlet boundary conditions specified on the entire boundary 
of the domain. Various versions of maximum principles and comparison principles are discussed 
in both continuous and discrete settings. In the literature, it is well-known that an acute-angled 
triangle is sufficient to satisfy the discrete weak maximum principle for pure isotropic diffusion. 
Herein, we show that this condition can be either too restrictive or not sufficient to satisfy various 
discrete principles when one considers anisotropic diffusivity, advection velocity field, or linear re¬ 
action coefficient. Subsequently, we derive appropriate restrictions on the mesh for simplicial (e.g., 
T3 element) and non-simplicial (e.g., Q4 element) elements. Based on these conditions, an iterative 
algorithm is developed to construct simplicial meshes that preserves discrete maximum principles 
using existing open source mesh generators. Various numerical examples based on different types 
of triangulations are presented to show the pros and cons of placing restrictions on a computational 
mesh. We also quantify local and global mass conservation errors using representative numerical 
examples, and illustrate the performance of metric-based meshes with respect to mass conservation. 


1. INTRODUCTION, MOTIVATION, AND CONTEMPORARY 

ADVANCEMENTS 

Diffusion-type equations are commonly encountered in various branches of engineering, sciences, 
and even in economics [1-3]. These equations have been well-studied in Applied Mathematics, and 
several properties and a priori estimates have been derived [4]. Numerous numerical formulations 
have been proposed and their performance has been analyzed both theoretically and numerically 
[5]. Several sophisticated software packages, such as ABAQUS [6], ANSYS [7], COMSOL [8], and 
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MATLAB's PDE Toolbox [9], have been developed to solve these types of equations. Special solvers 
for solving the resulting discrete equations have also been proposed and studied adequately [10]. 

This paper is concerned with numerical solutions for anisotropic advection-diffusion-reaction 
equations. Despite the aforementioned advances, it should be noted that a numerical solution 
always loses some mathematical properties that the exact solution possesses. In particular, the 
aforementioned software packages and popular numerical formulations do not satisfy the so-called 
discrete comparison principles (DCPs), discrete maximum principles (DMPs), and the non-negative 
constraint (NC). For example, consider the numerical simulation for pure anisotropic diffusion equa¬ 
tion in an L-shaped domain with multiple holes using the commercial software package ABAQUS. 

Numerical simulations are performed based on various unstructured finite element meshes (see 
Figure 1, which is to the scale), using the following popular anisotropic diffusivity tensor from 
hydrogeological and subsurface flow literature [11]: 

D(x) = RDeigenR^ (1-1) 


where Deigen is a diagonal matrix comprised of the eigenvalues of D (x). The corresponding principal 
eigenvectors are the column entries in the orthogonal matrix R. The expressions for Deigen and R 
are assumed as follows: 


R = 


D 


eigen 


cos{9) — sin(0) \ 
sin(0) cos{9) ) 

dmax 0 \ 

0 dmin J 


( 1 . 2 ) 

(1.3) 


Herein, dmax and dmin correspond to the maximum and minimum eigenvalues. 9 corresponds to the 
angle of orientation of the eigenvector coordinate system. It should be noted that these eigenvalues 
have physical significance and are related to the transverse and longitudinal diffusivities in the 
eigenvector coordinate system. Diffusion process is simulated based on equations (2.1a)-(2.1b). 
We assume v(x) = 0, a(x) = 0, and /(x) = 0. The prescribed concentration on the sides of 
the L-shaped domain is equal to zero. Correspondingly, the concentration on the perimeter of the 
holes are set to be equal to one. The values of dmax, dminj and 9 for D(x) given by the equations 
(1.1)-(1.3) are assumed to be equal to 1000, 1, and 

Very fine triangular (where the total number of nodes and mesh elements are equal to 86326 
and 169453) and quadrilateral (where the total number of nodes and mesh elements are equal to 
91778 and 90625) meshes are used to perform ABAQUS numerical simulations. The concentration 
profile obtained is shown in Figure 2. In this figure, we have not shown the concentration contour 
using four-node quadrilateral mesh, as it is almost identical to that of the contour obtained by 
employing three-node triangular mesh. The white area within the L-shaped domain with multiple 
holes represents the region in which the numerical value of concentration is negative and also exceeds 
the maximum. To be precise, in the case of triangular mesh, 2.22% and 3.92% of the nodes have 
violated the non-negative and maximum constraints. Correspondingly, the minimum and maximum 
values for concentration obtained are -0.0238 and 1.0076. These are considerably far away from 
the possible values, which are between 0 and 1. Similarly, for quadrilateral mesh, these values are 
slightly lower. Quantitatively, these are around 2.17% and 3.75%. But the minimum and maximum 
values of concentration (-0.0287 and 1.0086) are slightly higher than that of the triangular mesh. 
Additionally, it is evident that more than 6% of the nodes have unphysical negative values for the 
concentration. There are nodes for which the concentration exceeded the maximum constraint. 
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Furthermore, from Figures 3 and 4, it is evident that these values do not decrease with mesh 
refinement. In general, there are three possible routes to overcome such limitations and satisfy 
DCFs, DMPs, and NC; which we shall describe below. 

1.1. Strategy - I: Mesh restrictions. The first approach is to place restrictions on the 
mesh to meet maximum principles and the non-negative constraint. For isotropic homogeneous 
diffusivity, Ciarlet and Raviart [ 12 ] have shown that numerical solutions based on the single-field 
Galerkin finite element formulation, in general, does not converge uniformly. It should however be 
noted that the single-field Galerkin formulation is a converging scheme. Ciarlet and Raviart have 
also shown that a sufficient condition for single-field Galerkin formulation to converge uniformly for 
isotropic diffusion is to employ a well-centered three-node triangular element mesh with low-order 
interpolation. 

The obvious advantage of this approach is that one can use the single-field Galerkin formulation 
without any modification. The drawback is that an appropriate computational mesh may not exist 
because of the required restrictions on the shape and size of the finite element. For example, it 
is not an easy task (sometimes it not possible) to generate a well-centered triangular mesh for 
any given two-dimensional domain [13]. Note that requiring a mesh to be well-centered is a more 
stringent than requiring the mesh to be Delaunay. In fact, a well-centered mesh is Delaunay but 
the converse need not be true. 

In scientific literature, there are numerous commercial and non-commercial mesh generators 
that produce premium quality structured and unstructured meshes for various complicated do¬ 
mains. For instance, the survey paper by Owen [14] accounts for more than 70 unstructured mesh 
generation software products. But it should be emphasized that Owen [14] rarely mentions about 
non-obtuse, acute, and anisotropic A4-uniform mesh generators. However, it is evident from the 
above discussion that these types of meshes have a profound impact on solving various important 
physical problems related to diffusion-type equations. In recent years, there has been considerable 
effort in developing such types of mesh generators. For example, some open source meshing software 
packages which are relevant to mesh restrictions methodology are as follows: 

• Non-obtuse and acute triangulations in 2D: aCute [15-17] (a meshing software, which is 
based on Triangle [18]) 

• Anisotropic M-uniform triangulations in 2D: BAMG [19] in FreeFem-|— |- [20,21], BL2D [22] 

• Anisotropic A4-uniform triangulations in 3D: MmgSd [23] 

• Locally uniform anisotropic Delaunay meshes (surface, 2D, and 3D): CGALmesh [24-27] 

However, the use of these mesh generators in the area of numerical analysis and engineering, in par¬ 
ticular, to construct mesh restrictions for diffusion-type equations to satisfy DGPs, DMPs, and NG 
is hardly known. Recently, Huang and co-workers [28-30] used BAMG to generate anisotropic sim- 
plicial meshes to satisfy various discrete properties for linear advection-diffusion-reaction equations. 
But in their research works, the computational domains under consideration are uncomplicated. 

1.2. Strategy - II: Non-negativity, monotone, and monotonicity preserving for- 
mnlations. The second approach is mainly concerned with developing new innovative numerical 
methodologies based on certain physical and variational principles, so that they satisfy DGPs, 
DMPs, and NC. Broadly, these methods can be classified into the following three categories: 

• Non-negative formulations: A numerical method belongs to the class of non-negative 
formulations if the resulting numerical solution satisfies certain DMPs and NC. 
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• Monotone formulations: A numerical method is said to be monotone if the resulting 
numerical solution satisfies certain DMPs, DCPs, and NC. 

• Monotonicity preserving formulations: A numerical formulation is said to be monotonicity 
preserving if the resulting numerical solution does not exhibit spurious oscillations within 
itself. 

It needs to be emphasized that a non-negative formulation need not satisfy monotone conditions, a 
monotone numerical method need not be monotonicity preserving, and vice-versa. It is still an open 
research problem to develop a numerical formulation that meets all the aforementioned properties. 

Some popular formulations and notable research works in this direction are hnite difference 
schemes (FDS) [31,32], mimetic hnite difference methods (MFDM) [33,34], hnite volume meth¬ 
ods (FVM) [35-37], and hnite element methods (FEM) [38-41]. It should be noted that most 
of these techniques are inherently non-linear. For example, the optimization-based hnite element 
methodologies proposed by Nakshatrala and co-workers [11,41-43] enforce the desirable properties 
as explicit constraints. This is achieved by constructing variationally consistent constrained mini¬ 
mization problems for various numerical formulations. But one should note that this comes with 
an additional computational cost. However, Nakshatrala et al. [11] have shown through numerical 
experiments that the additional computational cost is less than 10% . 

1.3. Strategy - III: Post-processing methods. The third approach is about post-processing 
(PP) based methods. In literature, there are various types of PP methods, which can be used to 
recover certain discrete properties for diffusion-type equations. Some of the notable research works 
in this direction include: 

• Local and global remapping/repair methods [44,45] 

• Constrained monotonic regression based methods [46] 

• Cutoff methods (also known as the clipping methods) [47,48] 

• A combination of remapping/repair methods and cutoff methods [49,50] 

We shall now give a brief description of the pros and cons of these methods. Nevertheless, it is 
very difficult to apply these techniques to recover DCPs, DMPs, and NC for higher-order FEM 
methods, as the shape functions can change their sign within the element. In addition, most of the 
above methods do not have a variational basis. 

The remapping/repair techniques proposed by Shashkov and co-workers are designed to improve 
the quality of the numerical solutions, so that they satisfy some discrete properties. Even though 
these are efficient, conservative, linearity and bound preserving interpolation algorithms, it should 
be emphasized that they are mesh dependent. In addition, application of such algorithms to 
anisotropic diffusion equations to satisfy DMPs and NC are seldom [49,50]. 

The post-processing procedure based on a constrained monotonic regression problem proposed 
by Burdakov et al. [46] is locally conservative, bound preserving, and monotonicity recovering 
method. This is a constrained optimization-based PP method, wherein one needs to specify various 
constraints in order to fulfill certain discrete properties. It is applicable to FDS, FVM, and FEM. 
But in the case of EEM, this PP method is valid only for linear and multi-linear shape functions. In 
order to construct appropriate constraints for the optimization problem, one needs to know a prior 
information on the lower bounds, upper bounds, and monotonicity of the numerical solution for 
the physical problem. In general, obtaining the qualitative and quantitative nature of the solution 
to a given physical problem is not always possible. If such information on the monotonicity, lower 
and upper bounds for the numerical solution is not available, then this methodology reduces to the 
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standard clipping procedure. In addition, one should note that it is not always possible to satisfy 
DCPs using this constrained monotonic regression algorithm. For instance, a counterexample 
similar to the research work by Nakshatrala et al. (see [11, Section 4]) can be constructed to show 
that it does not satisfy a DCP. 

Finally, we would like to emphasize that a posterior cutoff analysis is a variational crime. In 
general, this method is neither conservative nor satisfies DMPs and DCPs. The primary objective 
of this method is to cut off the values of a numerical solution if it is less than a given number 
(which is the cutoff value). Hence, it is called the cutoff method. In the case of highly anisotropic 
diffusion problems and for distorted meshes, this method predicts erroneous numerical results 
[11,42]. By specifying the cutoff value to be zero, it is always guaranteed to satisfy NC through 
this methodology. In addition, if the nature of the solution is known a prior, then one can also 
prevent undershooting and overshooting of the numerical solution by chopping off those values. 

1.4. Main contributions and an outline of this paper. Herein, we shall focus on the hrst 
approach of placing restrictions on the mesh to meet desired mathematical properties. We shall 
derive sufficient conditions for restrictions on the three-node triangular and four-node quadrilateral 
finite elements to meet comparison principles, maximum principles, and the non-negative constraint 
in the case of heterogeneous anisotropic advection-diffusion-reaction (ADR) equations. The notable 
contributions of this paper are as follows: 

(i) We provide an in-depth review of various versions of comparison principles, maximum princi¬ 
ples, and the non-negative constraint in the continuous setting. 

(ii) We derive necessary and sufficient conditions on the coefficient (i.e., the “stiffness”) matrix 
to satisfy discrete weak and strong comparison principles. 

(iii) A relationship between various discrete principles within the context of mesh restrictions, 
numerical formulations, and post-processing methods is presented. 

(iv) We propose an iterative method to generate simplicial meshes that satisfy discrete properties 
using open source mesh generators such as BAMG, FreeFem-|--|-, and Gmsh. 

(v) Different types of non-dimensional quantities are proposed for anisotropic diffusivity, which 
are variants of the standard Peclet and Damkohler numbers. These quantities are extremely 
useful in numerical simulations and have not been discussed in the literature. 

(vi) Lastly, several realistic numerical examples are presented to corroborate the theoretical find¬ 
ings as well as to show the importance of preserving discrete principles. 

The remainder of this paper is organized as follows. In Section 2, we present the governing 
equations for a general linear second-order elliptic equation and discuss associated mathematical 
principles: comparison principles, maximum principles, and the non-negative constraint. Section 3 
provides several important remarks on the continuous and discrete properties of elliptic equations. 
In Section 4, we shall derive mesh restrictions for the three-node triangular element and the rect¬ 
angular element to meet the discrete versions of maximum principles, comparison principles, and 
the non-negative constraint. Finally, conclusions are drawn in Section 5. 

We will denote scalars by lower case English alphabet or lower case Greek alphabet (e.g., 
concentration c and density p). We will make a distinction between vectors in the continuum and 
finite element settings. Similarly, a distinction will be made between second-order tensors in the 
continuum setting versus matrices in the discrete setting. The continuum vectors are denoted by 
lower case boldface normal letters, and the second-order tensors will be denoted using upper case 
boldface normal letters (e.g., vector x and second-order tensor D). In the finite element context, 
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we shall denote the vectors using lower case boldface italic letters, and the matrices are denoted 
using upper case boldface italic letters (e.g., vector f and matrix K). Other notational conventions 
adopted in this paper are introduced as needed. 

2. LINEAR SECOND-ORDER ELLIPTIC EQUATION AND ASSOCIATED 

MATHEMATICAL PRINCIPLES 

Let n C be a open bounded domain, where “nd” denotes the number of spatial dimen¬ 
sions. The boundary of the domain is denoted by dil, which is assumed to be piecewise smooth. 
Mathematically, dil := Q — Q, where a superposed bar denotes the set closure. A spatial point 
is denoted by x £ fl. The gradient and divergence operators with respect to x are, respectively, 
denoted by grad[«] and div[«]. Let c(x) denote the concentration field. We assume that Dirichlet 
boundary condition is prescribed (i.e., the concentration is prescribed) on the entire boundary.The 
remainder of this paper deals with the following boundary value problem, which is written in terms 
of a general linear second-order differential operator in divergence form: 

C[c] := —div [D(x)grad[c(x)]] -|- v(x) • grad[c(x)] -|- a(x)c(x) = /(x) in (2.1a) 

c(x) = cP(x) on (912 (2.1b) 

where C denotes the second-order linear differential operator, /(x) is the prescribed volumetric 
source, a(x) is the linear reaction coefficient, v(x) is the velocity vector field, D(x) is the anisotropic 
diffusivity tensor, and cP(x) is the prescribed concentration. Physics of the problem demands that 
the diffusivity tensor (which is a second-order tensor) be symmetric. That is, 

D^(x) = D(x) Vx £ 12 (2.2) 

Remark 2.1. In mathematical analysis, the divergence form is a suitable setting for energy 
methods. However, some studies on maximum principles do employ the non-divergence form, which 
can be written as follows: 

nd r\2 o 

= E (2.3) 

2,JI = 1 2=1 

where the coefficient {P)ij, (q)i, and r(x), can be related to the physical quantities such as the 
diffusivity tensor, velocity field, and linear reaction coefficient. It should be, however, noted that 
the non-divergence form exists irrespective of differentiability of the diffusivity tensor. If D(x) is 
continuously differentiable, then there exists a one-to-one correspondence between the divergence 
form and the non-divergence form. In such cases, the operator C in the divergence form given by 
equation (2.1a) can be put into the following non-divergence form [51, Chapter 6]: 

C[c] = —D(x) • grad [grad[c(x)]] -|- (v(x) — div [D(x)]) • grad[c(x)] -|- a(x)c(x) (2.4) 

where we have used the following identity in combination with equation (2.2) to obtain equation 
(2.4) 

div [D'^(x)grad[c(x)] = D(x) • grad [grad[c(x)]] -|- div [D(x)] • grad[c(x)] (2.5) 

Based on the nature of the coefficients and connectedness of the physical domain, different 
versions of maximum and comparison principles exist in the mathematical literature [51-53]. As 
stated earlier in this paper, we shall restrict ourselves to the boundary value problem given by the 
equations (2.1a)-(2.1b). Further analysis pertaining to Neumann boundary conditions and mixed 
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boundary conditions within the context of maximum principles, comparison principles, and the non¬ 
negative constraint is beyond the scope of this paper, and one can refer to references [4,53-55]. 
We shall say that the operator C, is elliptic at a point x G if 

0< A„,in(x)^-^<^-D(x)^< A„,ax(x)^-^ (2.6) 

where Amin(x) and Amax(x) are, respectively, the minimum and maximum eigenvalues of D(x). The 
operator L is said to be strictly elliptic if there exists a constant Aq, such that 

0 < Ao < Amin(x) V X G n 

and uniformly elliptic if 

0<^Hl^<+oo VxGfl (2.8) 

In the studies on maximum principles, it is common to impose the following restrictions on the 
velocity field v(x) and the linear reaction coefficient a(x): 

a(x) >0 V X G n (2.9a) 

a(x) — Idiv [v(x)] >0 V x G 12 (2.9b) 

Q ^ ^ < Po < +00 V X G 11 and V i = 1, • • • ,nd (2.9c) 

where /3o is a bounded non-negative constant. If (D)jj and (v), are continuous in H, then the 
operator C is uniformly elliptic for any bounded subdomain CC H (which means that 11^ is 
compactly embedded in H) and the condition given in equation (2.9c) holds. The restrictions given 
in equation (2.9b) can be relaxed in some situations (e.g., see references [29,30]). But the constraint 
on a(x) given by equation (2.9a) cannot be relaxed. If a(x) < 0, then equation (2.1a) is referred 
to as an Helmholtz-type equation, which does not possess a maximum principle. From the theory 
of partial differential equations, it is well-known that the aforementioned boundary value problem 
given by equations (2.1a)-(2.1b) satishes the so-called (weak and strong) comparison principles, 
(weak and strong) maximum principles, and the non-negative constraint. For future reference 
and for completeness, we shall briefly outline the main results. For a more detailed mathematical 
treatment, one could consult references [4,51,52]. 


(2.7) 


Theorem 2.2 (Continuous weak and strict weak maximum principles). Let C he a uni¬ 
formly elliptic operator satisfying the conditions given by equations (2.9a)-(2.9c). In addition, let 
D (x) he continuously differentiable. Suppose that c(x) G C‘^{Ll) n C°(I1) satisfies the differential 
inequality £[c] < 0 m H, then the maximum of c{x) in H is obtained on dQ. That is, c(x) possesses 
the weak maximum principle (wMP ), which can be written as follows: 


max [c(x)] < max 0, max [c(x)] 
xgH L 


Moreover, if a{x) = 0, then we have the strict weak maximum principle (WMF ): 


( 2 . 10 ) 


m^ [c(x)] = max [c(x)] (2-11) 

xeO xG(9f1 


Proof. For a proof, see references [51,52]. 


□ 
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Theorem 2.3 (Continuous strong and strict strong maximum principles). Let the do¬ 
main Ll be simply connected. Given that c(x) satisfies wMP and the conditions given in Theorem 
2.2, then c(x) cannot attain an interior non-negative maximum in Ll unless it is a constant. This 
means that, c(x) possesses the strong maximum principle (sMP ) if the following hold: 

max [c(x)] = [c(x)] = m > 0 ^ c(x) = m in 0 (2.12) 

xeH 

Moreover, if a(x) = 0 and c(x) satisfies WMP, then we have the strict strong maximum principle 
(^SMP ) given as follows: 

max [c(x)] = m^ [c(x)] = m c(x) = m in 12 (2.13) 

xeH 

Proof. For a proof, see references [51,52]. □ 


Theorem 2.4 (Continuous weak and strong comparison principles). Let ci(x), C 2 (x) G 
(7^(12) n (7^(12). Suppose C he a uniformly elliptic operator satisfying the conditions given by the 
equations (2.9a)-(2.9c). Then L is said to possess 

• the weak comparison principle (wCP) i/ci(x) and C 2 (x) satisfies wMP, C\ci\ < C[c 2 \ in 
12, and ci(x) < C 2 (x) on 912, then the following holds: 

ci(x) < C 2 (x) Vx G 12 (2.14) 

• the strong comparison principle (sCP^ z/ci(x) and C 2 (x) satisfies sMP, C[ci] < C[c 2 \ in 
12, and ci(x) < C 2 (x) on 912, then the following holds: 


ci(x) < C 2 (x) Vx G 12 
Proof. For proof, see reference [52]. 


(2.15) 

□ 


Numerical formulations based on the finite element method, finite volume method, and finite 
difference method exist to solve the boundary value problem given by the equations (2.1a)-(2.1b). 
It is well-known that the framework offered by the finite element method is particularly attractive 
in obtaining accurate numerical results for elliptic partial differential equations. In particular, the 
single-field Galerkin formulation is a very popular finite element formulation. In this paper, we shall 
use the single-field Galerkin formulation to derive mesh restrictions. It should be, however, noted 
that restrictions imposed on a mesh may alter if an alternate numerical formulation is employed. 
But the overall procedure presented in this paper can be employed to derive mesh restrictions for 
other numerical formulations. 


2.1. Single-field Galerkin formulation. Let us define the following function spaces: 

C := {c(x) G I c(x) = cP(x) on 912} (2.16a) 

W := {tc(x) G I rc(x) = 0 on 912} (2.16b) 

where 77^(12) is a standard Sobolev space [51]. For weak solutions, the regularity of the diffusivity 
tensor can be relaxed as follows: 

J tr [D(x)'^D(x)] dl2 < -boo 


( 2 . 17 ) 


where tr[«] is the standard trace operator used in tensor algebra and continuum mechanics [56]. 
Given two helds a(x) and 6(x) on a set T), the standard L 2 inner-product over T> will be denoted 
as follows: 

{a]b)^ = j a(x) • 6(x) dV (2-18) 

V 

The subscript on the inner-product will be dropped if P = P. The single-held Galerkin formulation 
for the boundary value problem given by equations (2.1a)-(2.1b) can be written as follows: Find 
c(x) £ C, such that we have 

B{w]c) = L{w) V r(;(x) £ >V (2-19) 

where the bilinear form and the linear functional are, respectively, dehned as follows: 

B{w]c) := (tc;a(x)c) -|- (rc; v(x) • grad[c]) -|- (grad[tc];D(x)grad[c]) (2.20a) 

L{w) := {w,f{x)) (2.20b) 


2.2. Discrete single-field Galerkin formulation. Let the computational domain P be de¬ 
composed into “A^e^e” non-overlapping open sub-domains, which in the hnite element context will 
be elements. That is, 

Nele 

n = ij o'" (2.21) 

e=l 

The boundary of is denoted as := — P®. Let P^(n'^) denote the vector space spanned by 

linear polynomials on the sub-domain 12®. We shall dehne the following hnite dimensional subsets 


of C and W: 

:= |c'"(x) £C I c^ix) £ C°(12); c'‘(x)|^, £ P^(12®); e = !,■■■ ,Nele^ (2.22a) 

:= {u;^(x) £ W I w^{x) £ C°(f2); £ P^(L!®); e = !,■■■ ,iVek| (2.22b) 

A corresponding hnite element formulation can be written as follows: Find c^(x) £ C^, such that 
we have 

B{w^-,c^) = L{w’^) Vu;^(x)£W'^ (2.23) 

where B{w^;c^) and L{w^) are, respectively, given as follows: 

B{w^]c^) := ^rc^;a(x)c^^ -|- v(x) • grad[c^]^ -|- ^grad[u;^]; D(x)grad[c^]^ (2.24a) 

L{w^)-.= [w^-f{x)) (2.24b) 


Let “ui” denote the total number of degrees-of-freedom, “n/” denote the free degrees-of- 
freedom, and “up” be the prescribed degrees-of-freedom for the concentration vector. Obviously, 
we have nt = rif + rip. We assume that nt,np > 2. After hnite element discretization, the discrete 
equations for the boundary value problem take the following form: 


Kc = r 


(2.25) 


where K = [Kff\Kfp] is the stiffness matrix, c = 


1 T 




is the vector containing nodal 


concentration, and r = [rj]^ is the corresponding nodal volumetric source vector. The stiffness 
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matrices K, Kff, and Kfp are, respectively, of size n/ xnt, Uf xnf, and n/ x Up. Correspondingly, 
the nodal concentration vectors c, Cf, and Cp are of sizes nt x 1, UfXl, and UpXl. Similar inference 
is applicable to the load vector r. 

Before we state a discrete version of (weak and strong) maximum and comparison principles, we 
introduce the required notation. The symbols ^ and ^ shall denote component-wise inequalities 
for vectors and matrices. That is, given two (hnite dimensional) vectors a and b 

a <b means that ai <bi\/ i (2.26) 

Correspondingly, given two matrices A and B 

A < B means that (^)ij < ^ (2.27) 

Similarly, one can define the symbol and In the remainder of this paper, we will be 

frequently using the symbols 0 and O, which, respectively, denote a zero vector and a zero matrix. 

We shall now briefly outline the main results corresponding to the discrete weak and strong 
maximum principles in the form of dehnitions and theorems. Using these results, we shall discuss 
in detail about discrete comparison principles. However, it should be noted that Theorem 2.8 and 
its proof are new and have not been discussed elsewhere. We shall also present the necessary and 
sufficient conditions on the stiffness matrices K and Kfp to satisfy different versions of discrete 
maximum principles and the non-negative constraint. For more details, see references [57,58]. 

Definition 2.5 (Discrete maximum principles [58]). A numerical formulation is said to 
possess 

• the discrete weak maximum principle (DwMP ) if 

r ^ 0 implies max [c] < max [0, max [cp]] (2.28) 

• the discrete strict weak maximum principle (DWMP ) if 

r ^ 0 implies max [c] = max [cp] (2.29) 

• the discrete strong maximum principle (DsMP ) if it possesses DwMP and satisfies the 
following condition: 

r ^ 0, and max [c] = max [cf] = m >0 implies c = ml (2.30) 

• the discrete strict strong maximum principle (DSMP ) if it possesses DWMP and satisfies 
the following condition: 

r ^ 0, and max [c] = max [cf] = m implies c = ml (2-31) 

where max [•] denotes the maximal element of a vector and the symbol 1 is the vector whose com¬ 
ponents are all equal to 1. 

Theorem 2.6 (Necessary and sufficient conditions to satisfy DMPs). The stiffness ma¬ 
trix K given by equation (2.25) is said to possess 

• the discrete weak maximum principle (DwMPx ) if and only if all of the following condi¬ 
tions are satisfied: 

{a)K~j hO (b) - KjjKfp hO (c) - KjjKfpl ^ 1 
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(2.32) 


• the discrete strict weak maximum principle (TDWMP /^) if and only if all of the following 
conditions are satisfied: 

hO (b) - K-jKfp hO (c) - K~jKfpl = 1 (2.33) 

• the discrete strong maximum principle (DsMP ^^:) if and only if all of the following condi¬ 
tions are satisfied: 

{a)Kjj (b) - KjjKfp yO (c) - KjjKfpl X 1 or - KjjKfpl = 1 (2.34) 

• the discrete strict strong maximum principle (DSMP /^) if and only if all of the following 
conditions are satisfied: 

yO (b) - KjjKfp (c) - KjjKfpl = 1 (2.35) 

Proof. For a proof, see reference [58]. □ 

Definition 2.7 (Discrete weak and strong comparison principles). A numerical formu¬ 
lation is said to possess 

• the discrete weak comparison principle (DwCPj if it satisfies DwMP, and 

Cl :< C 2 on dll and ri :< r 2 in D implies Ci -< C 2 in D (2.36) 

• the discrete strong comparison principle (DsCP ) if it satisfies DsMP, and 

Cl :< C 2 on dH and ri -< r 2 in D implies ci -< C 2 in D (2.37) 

Theorem 2.8 (Necessary and sufficient conditions to satisfy DCPs). Let ci and C 2 
he two nodal concentration vectors corresponding to the volumetric source vectors ri and r 2 based 
on the equation (2.25). If ci and C 2 satisfy DwMP and the hypothesis o/DwCP (i.e., ci ^ C 2 
on dfl and ri ^ r 2 in LI), then a necessary and sufficient condition to satisfy the discrete weak 
comparison principle (DwCP^J (which means that ci < C 2 in LI) is that the stiffness matrix K 
possess DwMP/<: (which is given by equation (2.32) in Theorem 2.6). 

If Cl and C 2 satisfy DsMP and the hypothesis o/DsCP, (i.e., ci ^ C 2 on dLl and ri -< r 2 in 
LI), then a necessary and sufficient condition to satisfy the discrete strong comparison principle 
(DsCPjcj (which means that ci ^ C 2 in Ll) is that the stiffness matrix K possess DsMPx (which 
is given by equation (2.34) in Theorem 2.6). 

Proof. For convenience, let us define the following: 

C 3 := Cl - C 2 (2.38a) 

7*3 := ri - r2 (2.38b) 


Clearly, C 3 and r 3 satisfy the following: 


Kc 3 = r3 


II 


(2.39) 


Necessary condition to satisfy DwCP/c : Let ci ^ C2 in O, which implies that C3 ^ 0 in Q. The 
hypothesis of DwCPjc and the fact that C 3 ^ 0 in imply the following: 


T*3 ^ 0 

in H 


(2.40a) 

C 3 ^ 0 

on (9H 


(2.40b) 

max [C 3 ] <0 

on dQ 


(2.40c) 

max [C 3 ] < max 
n 

0, inax [C 3 ] 
oil 

= 0 

(2.40d) 


which implies that C 3 satisfies DwMP (based on equation (2.28) in Definition 2.5). But vector C 3 also 
satisfies equation (2.39). Hence, according to equation (2.32) and the hypothesis of Theorem 2.6, 
it is evident that K must possess DwMP/^. This completes the proof for the necessary condition 
to satisfy DwCPjc- 

Sufficient condition to satisfy DwCPx : It is given that ci and C 2 satisfy DwMP. Equations 
(2.38a)-(2.38b) and DwCPjc imply that 


C 3 ^ 0 on (9H (2.41a) 

r 3 ^ 0 in H (2.41b) 


If the stiffness matrix K possess DwMP^, it is evident from Theorem 2.6 and equations (2.39), 
(2.41a)-(2.41b) that vector C 3 satisfies DwMP. Hence, from Definition 2.5 and equations (2.28), 
(2.41a)-(2.41b), we have the following result: 


nmx [C 3 ] < max 

n 


0, inax [C 3 ] 
oil 


= 0 


(2.42) 


From equation (2.42), it is evident that the least upper bound for any component of vector C 3 is 
equal to zero. Hence, we have C3 ^ 0 in H. This implies that ci ^ C2 on H, which completes the 
proof for the sufficient condition to satisfy DwCPk- 

Necessary condition to satisfy DsCPk : Following the arguments about the proof for the nec¬ 
essary condition to satisfy the DwCPk property, it is evident that C 3 satisfies DwMP. In addition, 
we are given that ci C 2 in H. This implies C 3 -< 0 in H. Based on the hypothesis of DsCPk and 
utilizing the fact that C 3 -< 0 in H yields the following: 


max [C 3 ] < 0 on (2.43a) 

max [C 3 ] <0 in H (2.43b) 

This means that vector 0 is the least upper bound for C 3 in H, and any component of C 3 is strictly 
less than zero in the interior of the domain H. From equation (2.43a) and (2.43b), it is clear that 
the non-negative maximum value for vector C 3 occurs on the boundary dQ. From equation (2.30) in 
Definition 2.5, it follows that vector C 3 satisfies DsMP. Hence, according to conditions specified by 
equation (2.34) and the hypothesis of Theorem 2.6, it is evident that K must possess the DsMP^ 
property. This completes the proof for the necessary condition to satisfy DsCP/^. 

Sufficient condition to satisfy DsCPk : Given that ci and C 2 satisfy DsMP. Under the assump¬ 
tions of DsCPic and from equations (2.38a)-(2.38b), we have the following relations: 

rs -< 0 in H 

C 3 ^ 0 on (9H 
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(2.44a) 

(2.44b) 






If the stiffness matrix K possess DsMP/^, it is evident from Theorem 2.6 and equations (2.39), 
(2.44a)-(2.44b) that vector C 3 satisfies DsMP. Hence, by appealing to Definition 2.5 and equation 
(2.30), if vector C 3 does not attain a non-negative maximum value at an interior point of D, then 
we have the following result: 


max [C 3 ] < max 


0, inax [C 3 ] 
oil 


= 0 


(2.45) 


which implies that each component of vector C 3 is less than zero. Hence, we have ci C 2 in 
D. Suppose, if vector C 3 attains a non-negative maximum value at an interior point of D, then 
according to the surmise of DsMP, we first need to satisfy DwMP. So from equation (2.28), we 
have the following relation: 


rnnx [ 03 ] < max 

o 


0 , max [C 3 I 

an 


= 0 


(2.46) 


Secondly, according to DsMP, we also need to satisfy the equation (2.30). These conditions in 
terms of vector C 3 are given as follows: 


nmx [cs] = max [03] = m > 0 (2.47a) 

C 3 = ml in Q (2.47b) 

From equations (2.46) and (2.47a)-(2.47b), it is evident that m = 0; which implies that C3 = 0. 
Thus, we have ci = C 2 in D. But from equation (2.39), it is obvious that rs = 0 in D, which 
contradicts the hypothesis of DsCPk given by the equation (2.44a). Hence, we have the final result 
Cl ^ C 2 in D, which completes the proof for the sufficient condition to satisfy DsCPk. □ 


In the next section, we shall discuss the various factors (i.e., mesh restrictions, numerical 
formulations, and post-processing methods) that influence the satisfaction of discrete versions of 
maximum principles, comparison principles, and the non-negative constraint. 


3. YET ANOTHER LOOK AT CONTINUOUS AND DISCRETE PRINCIPLES 

Based on the finite element methodology outlined in subsection 2.2, we shall analyze the prop¬ 
erties that the stiffness matrix K inherits from the continuous problem. An important attribute 
that the discrete system needs to have in order to mimic the mathematical properties that the 
continuous system possesses is that the stiffness matrix K// has to be a (reducible or irreducible) 
monotone matrix. The part (a) in all the equations (2.32)~-(2.35) of Theorem 2.6 corresponds to 
reducibility or irreducibility of K//- 

On general computational grids, it is well-known that the stiffness matrix K// obtained via 
low-order finite element discretization might not be a monotone matrix [12, 39, 58]. So, the 
discrete single-field Galerkin formulation might (or shall) violate the non-negative constraint, 
discrete maximum principles, and discrete comparison principles on unstructured computational 
meshes [12,40-42,57,58]. The violation is more severe if the diffusion tensor is anisotropic. One 
of the ways to overcome such unphysical values for concentration and preserve the discrete prop¬ 
erties is to restrict the element shape and size in a computational mesh. This can be achieved by 
developing sufficient mesh conditions under which K fj is ensured to be a reducibly or irreducibly 
diagonally dominant matrix [59]. Before we discuss such a class of monotone matrices, which are 
easily amenable for deriving mesh restrictions, some important remarks on various DMPs and their 
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relationship to DCPs and NC are in order. We would like to emphasize that such a comprehensive 
discussion is not reported elsewhere in the literature. 

3.1. Simply connected vs. multiple connected domains. For many applications in 
mathematics, sciences, and engineering, it is necessary to at least satisfy the weak or strict weak 
maximum principle. But there are numerous cases where in it is required to satisfy a strong version 
of the maximum principle [39,57]. In such scenarios, geometry and topology of the domain play 
a vital role. According to the hypothesis of Theorem 2.3, it is evident that a strong maximum 
principle exists if the domain is simply connected (see reference [52, Chapter 3]). However, one 
should not immediately conclude that if a domain is not simply connected, then a strong maximum 
principle will not exist [53,58]. 

In a discrete setting, Ishihara [57], Draganescu et.al. [39], and Mincsovics and Hovarth [58] have 
conducted various numerical experiments related to discrete strong maximum principles for multiple 
connected domains. They performed analysis related to satisfaction of DsMP/^ and DSMP^ 
for various non-obtuse and acute triangulations for multiple connected domains. In particular, 
Mincsovics and Hovarth discuss various interesting examples related to the irreducibility property 
of the stiffness matrix Kff when the domain is not simply connected. In all of their examples, 
they solve the following equations: 

ac — Ac = 0 in H (3.1a) 

c(x) = cP(x) on do. (3.1b) 

where the linear decay a = 0 or a = 128. Through numerical experiments, the authors demon¬ 
strate that even though the triangulation satisfies the non-obtuse or acute angled mesh condition 
(proposed by Ciarlet and Raviart [12]), it is not guaranteed to fulfill either DsMP or DSMP. This 
means that non-obtuse [15] and well-centered triangulation [13] of any given domain will always 
satisfy the weak DMPs, but need not satisfy the strong DMPs. 

Within the context of directed graphs [39,59], there is a one-to-one correspondence between 
irreducibility of the stiffness matrix Kff and the interior vertices of the computational mesh [30]. 
In order to satisfy the discrete (strong and strictly strong) maximum principle, the mesh has to 
be interiorly connected, which in turn implies that Kff has to be irreducible [60]. By interiorly 
connected mesh, we mean that any pair of interior vertices of the mesh are connected at least by 
an interior edge path [39]. Hence, Kjj >- 0 and —K~^jKfp 0 in Theorem 2.6 correspond to 
this discrete connectedness property of the computational mesh [39,58]. However, it should be 
noted that irreducibility is a necessary condition, but not sufficient. For other details on numerical 
aspects related to mesh connectivity, see references [58, Section 4, Figures 1-4], [39], and [30]. 

3.2. Minimum principles, and non-negative and min-max constraints. Due to linear¬ 
ity of the operator C, similar theorems corresponding to minimum principles and the non-negative 
constraint for equations (2.1a)-(2.1b) can be derived. To obtain the non-negative solution and corre¬ 
sponding min-max constraint on c(x), we shall appeal to the continuous weak minimum/minimum- 
maximum principle, which can be written as follows [51,52]: 

Lemma 3.1 (Continuous weak minimum/minimum-maximum principle). Let C be a uniformly 
elliptic operator satisfying the conditions given by (2.9a)-(2.9c) and D(x) be continuously differen¬ 
tiable. Given that C[c] > 0, cP(x) > 0, and c(x) £ (7^(12) n C^{Q), then c(x) possess a continuous 
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weak minimum principle, which is given as follows: 


min [c(x)] > min 0, min [c(x)] 
xgQ [ xean 


(3.2) 


Moreover, if C[c\ = 0, then we obtain the classical weak minimum-maximum principle for c(x) in 
Q, which is given as follows: 


min 

0, min [c(x)] 

< c(x) < max 

0, max [c(x)] 




xg90 


(3.3) 


Proof. For a proof, see references [51,52]. 


□ 


It is evident from equation (3.2) that for /(x) > 0 and cP(x) > 0, we have c(x) > 0 for 
any x S 0. Correspondingly, a discrete version of continuous weak minimum principle and weak 
minimum-maximum principle is given as follows: 


Definition 3.2 (Discrete weak minimum/minimum-maximum principle). A numerical for¬ 
mulation is said to possess 

• the discrete weak minimum principle if 

r ^ 0 implies min [c] > min [0, min [cpj] (3.4) 

• the discrete weak minimum-maximum principle if 

r = 0 implies Cmml ^ C ^ Cmaxl (3.5) 

where Cmin := Riin[cp], Cmax := max[cp], and min [•] denotes the minimal element of a vector. 

3.3. High-order finite element methods. An attribute of low-order finite elements, which 
plays a central role in designing non-negative formulations for diffusion-type equations, is that the 
shape functions for these elements are monotonic and do not change their sign within the element 
[42]. Moreover, they are convenient to generate computational meshes for complex geometries [61], 
to perform error analysis [62], and for adaptive local mesh refinement [63]. High-order hnite 
elements are widely use for solving smooth problems, as one can obtain exponential convergence 
under high-order interpolations for these problems. But the shape functions of high-order finite 
elements change the sign within an element, which makes them not suitable under most of the 
current non-negative formulations (e.g., see reference [42,64]). The conditions presented in this 
paper will also not be applicable to high-order finite elements for the same reason of change in sign 
of interpolation functions within an element. 

As compared to low-order finite element methods, the discrete counterparts of continuous weak 
and strong maximum principles for high-order finite element methods is not well understood yet. 
This is because of the complicated task related to the test of non-negativity of a multivariate 
polynomial [65]. It should be noted that construction of non-negative high-order shape functions 
for finite element methods is still an unsolved problem and its roots can be traced back to the 
famous Hilbert’s problem [66,67]. Within the context of variational methods, probably, the 
research works by Ciarlet [68,69] are the hrst attempt to develop high-order non-negative shape 
functions to satisfy a discrete maximum principle. This study is based on a general theory of 
discrete Green’s function (DGF) for uniformly elliptic linear partial differential operators. Later, 
various attempts were made by different researchers to develop shape functions and derive mesh 
restrictions based on the DGF approach. Most of them are pertinent to one-dimensional problems 
or particular cases of isotropic diffusion. The conditions to be met for high-order elements to satisfy 
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DMPs based on DGF methodology are much more stringent and have a less broad scope for general 
applications. Furthermore, one should be aware that the discrete analogues of continuous Green’s 
functions are applicable only for linear problems, and cannot be extended to to non-linear problems 
(such as semi-linear and quasi-linear elliptic partial differential equations). Hence, such a method 
will have limited scope. For more details, one can consult the following references [65,70,71]. 

3.4. Relationship between various DCPs and DMPs. It is evident from Dehnitions 2.5, 
2.7, and 3.2; Lemma 3.1; and Theorems 2.6 and 2.8 that if a numerical formulation satisfies either 
DwCP/DwCPic or DsCP/DsCPjC) then it automatically obeys DwMP/DwMP^ and NC. Figure 5 
illustrates a graphical representation among various numerical solution spaces that satisfy different 
DMPs and DCPs within the context of mesh restrictions. By a (finite dimensional) numerical 
solution space Vp, we mean a set of numerical solutions {cj}”^p, which satisfy a given discrete 
property given by p. For example, if a concentration vector Cj corresponding to a given volumetric 
source vector satishes the discrete property DwMP, then c* S VdwMP- It should be noted that 
within the context of DMPs, DCPs, and NC; Vp is not a vector space. This is because if c* £ Vp, 
then according to non-negative property —Cj ^ Vp, which is one of the properties needed for Vp to 
be a vector space. 

Herein, we would like to emphasize that the route taken to satisfy p is very important. One 
can fulhll a discrete property p in numerous ways. In general, this is achieved by either placing 
mesh restrictions or developing a new (non-negative or monotone or monotonicity based) numerical 
formulation or through various post-processing methods. Couple of these techniques are developed 
based along the lines similar to Theorem 2.6 and others based on Dehnition 2.5. But it should be 
noted that developing numerical formulations accordant to Theorem 2.6 is much more difficult than 
that of Definition 2.5. This is because in order to satisfy Theorem 2.6, we need to place restrictions 
on the stiffness matrices R// and R/p. On the other hand, the hypothesis of Definition 2.5 does 
not assume any particular constraints on K. Hence, we would like to differentiate between the set 
of discrete properties given by DwMP^, DWMP/^, DsMP^, DSMP^, DwCP/^, and DsCPk to 
that of DwMP, DWMP, DsMP, DSMP, DwCP, and DsCP. 

In spite of the fact that there are several numerical methods available to satisfy a given discrete 
property p, from the characterization of Vp, it is evident that the resulting numerical solution 
spaces will be the same (for example, we have VowCPk = I^DwCp)- From Theorems 2.6 and 2.8, it 
is evident that among various DMPs and DCPs, we have the following set inclusions: 


VdsmPjc C VdsMPjc C VowMPif C VdwMPjc (3.6a) 

VdsCP/c VdwCPjc (3.6b) 

But it should be noted that a similar type of enclosure for numerical solution spaces between DMPs 
and DCPs does not hold: 

VDwCPif ^ VdwmPjc (3.7a) 

VowCPif 't VdsMP/c (3.7b) 

VDwCPic 't VdSMPjc (3.7c) 

VosCPic ^ I^DSMPjf (3.7d) 


The reason for such a non-enclosure stems from the hypothesis of Theorems 2.6 and 2.8, wherein 
we only need to satisfy DwMP/^ for DwCP/^ and DsMP/^ for DsCP/^. In a discrete setting, a 
numerical methodology may inherit one or more than one of these discrete principles and in some 
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cases none. Now, we shall discnss in detail a class of numerical formulations, which satisfy a certain 
discrete property p. We shall epitomize our findings based on various popular research works in 
literature, which span across different disciplines, such as computational geometry, optimization 
theory, numerical linear algebra, and partial differential equations. 

(a) Isotropic and anisotropic non-obtuse angle conditions: Recently, Huang and co-workers [28-30] 
were able to satisfy certain discrete properties through mesh restrictions, which are based on 
anisotropic A^-uniform mesh generation techniques. However, their theoretical investigation is 
mainly restricted to linear simplicial elements, specifically, the three-node triangular element 
and four-node tetrahedral element. But one should note that it is difficult to extend the 
procedure outlined by Huang and co-workers to multi-linear elements, such as the four-node 
quadrilateral element, six-node wedge element, and eight-node brick element. This is because 
the partial derivatives of the shape functions for multi-linear finite elements are not constant 
(for more details, see subsection 4.4 of this paper, which discusses mesh restrictions for a 
rectangular element). Nevertheless, constructing a WCT mesh [13] or an anisotropic A4- 
uniform triangular mesh [72,73] that satisfies various DMPs for an arbitrary domain is still an 
open problem [74]. The key-concept we would like to emphasize is that the numerical solutions 
obtained for isotropic diffusion-type equations using WCT meshes and anisotropic diffusion- 
type equations nsing the diffnsivity tensor based anisotropic A4-uniform meshes satisfy all 
versions of discrete maximum principles. In addition, if the hypothesis of DwCPk and DsCP^ 
is satisfied, then these meshes also satisfy all versions of discrete comparison principles. 

(b) Non-linear finite volume and mimetic finite difference methods: Le Potier’s method [35] and 
Lipnikov et al. [75] are some of the noteworthy works in the direction of FVM that satisfy 
the non-negative constraint, but do not possess a discrete version of the comparison principles 
and the maximum principles. However, it shonld be noted that recently these authors have 
developed techniques based on non-linear finite volume methods [35, 37] and mimetic hnite 
difference methods [34] to satisfy various versions of DMPs for a certain specific class of linear 
self-adjoint elliptic operators. But it should be noted that there is no discnssion on satisfying 
various DCPs. 

(c) Optimization-based finite element methods: Based on the works by Liska and Shashkov [40] 
and Nakshatrala and co-workers [11,41-43], the optimization-based low-order finite element 
methods, under certain conditions (when Kfj is symmetric and positive definite), can be 
written as follows: 


^ff^f — ^ ^fP^P ^min Amax 

(3.8a) 

C* ■ 1 ^ Cf ^ C* 1 

'-'mm-*- — — '-'max-*- 

(3.8b) 

Amin ^ 0 

(3.8c) 

Amax ^ 0 

(3.8d) 

{Cf — Cminl) ■ Amin = 0 

(3.8e) 

(Cmaxl ~ ^/) ■ -^max = 0 

(3.8f) 


where and are the minimum and maximum concentration values possible in H. These 
values can be obtained based on the bonndary conditions and a prior knowledge about the 
solution. Amin is the vector of Lagrange multipliers corresponding to the constraint cN^l Cf 
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and similarly Amax is the vector of Lagrange mnltipliers corresponding to the constraint cj ^ 
c* 1 

^max * 

Based on the nature of constraints, one can satisfy different discrete principles and it should 
be emphasized that DsMP, DSMP, DwCP, and DsCP can be fulfilled only under certain con¬ 
ditions. If either r 0 or r 0, then the non-negative constraint and the weaker versions of 
discrete minimum/maximum principles can be satisfied by specifying either or But 

in the case of r = 0, both and can be prescribed based on the Dirichlet boundary con¬ 
ditions; moreover, according to Definition 3.2 and from equation (3.5), we have = Cmin and 
Cmax = Cmax- However, One should note that if the qualitative and quantitative nature of the 
solution is known a prior, then one can satisfy all of the discrete versions of (weak and strong) 
maximum principles by specifying and ia the Karush-Kuhn-Tucker conditions given 
by equations (3.8a)“(3.8f). In general, these methods do not inherit a discrete strong maximum 
principle and a discrete comparison principle. For more details, a counter example is shown 
in the reference [11, Section 4, Figure 1]). Nevertheless, satisfying DsMP, DSMP, DwCP, and 
DsCP is still an open problem and are interesting topics to investigate in future endeavors. 

(d) Variationally inconsistent methods: In literature, there are various post-processing methods 
[46-48] available that can recover certain discrete properties if a prior information about the 
numerical solution is known. However, one should note that such methods are variationally 
inconsistent. A summary of the above discussion between various discrete principles within the 
context of FDS, MFDM, FVM, FEM and PP based methods is pictorially described in Figure 
6 . 

In the next section, we shall derive sufficient conditions on the three-node triangular element 
and four-node quadrilateral element to satisfy discrete versions of comparison principles, maximum 
principles, and the non-negative constraint. 

4. MESH RESTRICTIONS TO SATISFY DISCRETE PRINCIPLES 

In this section, we shall utilize and build upon the research works of Huang and co-workers 
[28-30, 72] for linear second-order elliptic equations. We first present, without proofs, relevant 
mathematical and geometrical results required to obtain mesh restrictions for simplicial elements. 
We will then use these results to construct mesh restriction theorems and generate various types of 
triangulations (see Algorithm 1 and the discussion in subsection 4.3) using the open source mesh 
generators, such as Gmsh [76] and BAMG [19] available in the FreeFem-|— |- software package [20,21]. 
It should be emphasized that these results cannot be extended to Q4 element, as this element is not 
simplicial. For more details on mesh restrictions for Q4 element, see subsection 4.4 of this paper. 

Let a;i, X 2 , ■ ■ ■, Xnd+i denote the vertices of an arbitrary simplex Hg S Th, where 7h is a sim- 
plicial triangulation of the domain D. The subscript ‘/i’ in the triangulation Th corresponds to 
the maximum element size (which will be described later in this section, see equation (4.20a)). 
Based on various values of h, we have an affine family of such simplicial meshes denoted by {Th}- 
Designate the total number of vertices and the corresponding interior vertices of Th by ‘An’ and 
‘Niv\ The edge matrix of He, which is denoted by is dehned as follows: 

EQ^:=[x2-Xi,X3-Xi,...,Xnd+l-Xl] y^e^Th (4.1) 

then an edge connecting vertices Xp and Xq of Hg is denoted as epq. Correspondingly, the edge 
vector Bpq^fi^ (which can be expressed as a linear combination of the elements corresponding to the 
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edge matrix Efie) the element boundary dQe in-terms of Xp, Xg, and Cpg are given by: 

epq,fi, = Xq - Xp '^p,q = l,2,...,nd+l and p/g (4.2a) 

nd-\-l 

dn^ = U epq (4.2b) 

p,q=l 
p¥=g 

Following references [28,77], a set of q-vectors corresponding to this edge matrix are defined 
as follows: 


nd-\-l 

■ 1 Qnd-\-l\ 

(4.3a) 

(li+ Qp = 0 

p=2 

VDg G % 

(4.3b) 


Let (fp^ denote the linear basis function associated with the Pg-th. global vertex in the triangulation 
Th- Then, {(Ppg}p^=i and {ppg}pJ'Zi span the respective finite dimensional subsets and given 
by equations (2.22a)-(2.22b). Denote the face opposite to vertex Xp by Fp and the corresponding 
unit inward normal pointing towards the vertex Xp by rip. The perpendicular distance (or the 
height) from vertex Xp to face Fp is denoted by hp. In the case of 2D, the above set of geometrical 
properties are pictorially described in Figure 7. 

Definition 4.1 (Positive linear maps). Let Win ■= Mnxn(lK) be the set of all real matrices 
of size n X n, whieh forms a vector space over the field M. A linear map ‘h : —)• is called 

positive if^{A) is positive semi-definite whenever A is positive semi-definite. Similarly, is called 
strictly positive if ^{A) is positive definite whenever A is positive definite. 

Theorem 4.2 (Strictly positive linear mapping of anisotropic diffusivity). Tef4>[«] 
Jq [•jdD. Show that $[•] is a linear map and 4>[D(x)] is symmetric, uniformly elliptic, and bounded 
above. 


Proof. From Definition 4.1, it is evident that ‘h[«] is a linear map and <h[D(x)] is symmetric. 
From equations (2.6) and (2.21), we have ^ • D(x)^ > 0 and meas(De) > 0. It is well known that 
Lebesgue integration of a scalar for a strictly positive measure is always greater than zero. Hence, 
^ ■ D(x)^dD > 0 Vx G Dg. Now, integrating equation (2.6) over Dg results in the following 
relation: 

0 < 7min^ • ^ <- 7777 [ ^ ■ D(x)^dD < ■ $, G and Vx G Dg (4.4) 

meas(sZg) J 
fie 


where the positive constants 7min and 7max are respectively the integral average of minimum and 
maximum eigenvalues of D(x) over Dg. These are given as follows: 


7min • — 


meas(De) 


Amin(x)dD 


1 


Tmax • — 


meas(Dg 


Amax(x)dfl 


(4.5a) 


(4.5b) 
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Since the vector ^ is independent of x and Vie, we can interchange the order of integration. This 
gives us the following equation: 

0 < 7min^ • ^ ^ 7max4 G and Vx G Vie (4.6) 

niG3jS(^S 2g j 

which shows that d>[D(x)] is a strictly positive linear map of D(x) and indeed preserves its prop¬ 
erties. □ 


4.1. Geometrical properties and finite element analysis of simplicial elements. Based 
on the above notation, we have the following important mathematical results relating the linear 
basis functions, hnite element matrices, and geometrical properties of simplicial elements. 

(1) The integral element average anisotropic diffusivity is given by: 

fie 


is a strictly positive linear map (see Definition 4.1 and Theorem 4.2). 

(2) For any arbitrary simplicial element Vie V Th, the vector associated with the face Fp, 
the gradient of the linear basis function grad[(/?pg], the unit inward normal rip, and the 
height hp are related as follows [74,77]: 


Qp = grad[99pj 


^ Vp= 1,2,... ,rad-|- 1 

hp 


(4.8) 


(3) The dihedral angle (3pq (measured in the Euclidean metric) between any two faces Fp and 
Fq is related to q-vectors (qp and q^) and unit inward normals (rip and rig) by the following 
equation [28,29]: 


cos(/3pg) = -rip • rig = 


qp • qg 

l|qplll|qgll 


Vp, q = 1,2,... , nd + 1 and p/q 


(4.9) 


where 


is the standard Euclidean norm. Similarly, the dihedral angle /3 ^-i measured 


■-1 


in metric is given as follows: 


cos(/3 ~-i) = 
^ vqFne 


qp Dflefjq 


Vp, q = 1,2,..., nd + 1 and p ^ q 


(4.10) 


where || • denotes the norm in metric. For example, llqpll^^^ = y qp ' -^^Oeqp- 

(4) For any simplex Vie vTh, the gradient of the linear basis functions (grad[(ppg] and grad[(pgg]), 
the dihedral angle jSpq, and heights {hp and hg) are related as follows [28,29]: 


meas(De) grad[ppj • grad[(pgj 


meas(De) cos(/3j 


•pqj 


hphq 


Vp, q = 1,2,..., nd + 1 and p^q 


(4.11) 


and in 2D, equation (4.11) reduces to: 


meas(De) ( grad[p. 


-'PsJ 


grad[p, 


Qgi 


COt(/3pg) 


Vp, q = 1, 2,..., nd-|-1 and p^q 
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(4.12) 














(5) For any arbitrary simplicial element fig £ Th, the integral element average anisotropic 
diffusivity the gradient of the linear basis functions (grad[(/?pg] and grad[(/jqg]), the 

dihedral angle j3 measured in -D^e ™6tric, and heights {hp and hq) are related as 

follows [29]: 


meas(17e) ( grad[99pj • grad[(^qj 


and in 2D, equation (4.13) reduces to: 


meas(De) ( grad[(pp^] ■ grad[(^gj 


meas(ilp) cos(B 


1-1 


KIIb„ iiq, 


1-1 

Dn. 


\/p,q = 1,2,... ,nd+I and p^q (4-13) 


detfDt 


■cot(/3 ~-i) 

'\/p,q = 1,2,... ,nd + 1 and p/g (4-14) 


4.2. Sufficient conditions for a three-node triangular element. Using the mathematical 
results outlined in subsection 4.1, we shall present various sufficient conditions on the T3 element 
to satisfy different types of discrete properties. In general, there are two different approaches to 
obtain sufficient conditions. The first approach, which shall be called global stiffness restriction 
method, involves manipulating the entries of the global stiffness matrix, so that it is either weakly 
or strictly diagonally dominant (see Theorems 4.3 and 4.4) based on the nature of a(x). This 
means that Kfj exists and K 7 / - 0. The component wise entries of Kfj satisfy the following 
conditions: 

(a) Positive diagonal entries: {Kfj).. > 0, 

(b) Non-positive off-diagonal entries: <0 Vz / j, and 

one of the following two conditions: 

(c) Strict diagonal dominance of rows: | (Ffjj) .. | > 

(c) Weak diagonal dominance of rows: | (Kff).. \ > '''■J 

The second method, which shall be called local stiffness restriction method, engineers on stiffness 
matrices at the local level, so that they are weakly diagonally dominant. Once we ascertain that 
all of the local stiffness matrices are weakly diagonally dominant, then the standard finite element 
assembly process [78] guarantees that the global stiffness matrix Kff is monotone and weakly 

diagonally dominant. In particular, if a(x) > 0, then Kjf is strictly diagonally dominant. It 

should be noted that the above three conditions are sufficient, but not necessary, for the global 
stiffness matrix Kjf tohe monotone. 

4.2.1. Global and local stiffness restriction methods. In general, to get an explicit analyt¬ 
ical formula for Kjj is extremely difficult and not practically viable. Hence, it is not feasible to find 
mesh restrictions based on the condition that Kjj ^ 0. So an expedient route to obtain monotone 
stiffness matrices through mesh restrictions is by means of weakly or strictly diagonally dominant 
matrices, which form a subset to the class of monotone matrices [60, Section 3, Corollary 3.20 
and Corollary 3.21]. The obvious edge being that there is no need to compute {Kjj)ij explicitly. 
Based on the global stiffness restriction method, we shall now present stronger and weaker mesh 
restriction theorems. A brief and concise proof to the theorems shall be put forth, which is along 
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similar lines to the one presented in references [29,30]. These mesh restriction theorems shall be 
used in constructing triangular meshes to satisfy different discrete principles (see subsection 4.3). 

Theorem 4.3 (Anisotropic non-obtuse angle condition). If any nd-simplicial mesh sat¬ 
isfies the following anisotropic non-obtuse angle condition: 

hri 11 U 11 hr) hn 11 CK 11 7^ 

(nd + l)(nd + 2)A_ 5^^ - 

Vp, g = 1, 2, • • • , nd + 1, pfiq, He^Th (4.15) 


then we have the following three results: 


• The global stiffness matrix Kff is (reducibly/irreducibly) weakly diagonally dominant if 
a(x) > 0 and is (reducibly/irreducibly) strictly diagonally dominant if a{x) > 0 for all 

X G n. 

• The discrete single-field Galerkin formulation given by equations (2.24a)-(2.24b) in com¬ 
bination with equations (4.7)-(4.14) satisfies DwMP^/DWMPx; where and 

ll“llcx)T7e defined as: 

||u||cx 3 q ;= maximize ||u(a;)|| (4.16a) 

’ ^ X£0,e 

||a||ooQ := maximize a{x) (4.16b) 

’ cc€f!e 


and denotes the minimum eigenvalue of D^i^; hp, 


hn, and B ~-i are respectively 


the heights and metric based dihedral angle opposite to the face Fr of element hie- 
• Moreover, if the triangulation Th is interiorly connected, then the global stiffness matrix 
Kff is irreducibly weakly or strictly diagonally dominant based on the nature of a{x) and 
the diserete single-field Galerkin formulation given by equations (2.24a)-(2.24b) satisfies 
DsMPk/DSMPk. 


Proof. For proof, see References [29,30]. 


□ 


Theorem 4.4 (Generalized Delaunay-type angle condition). In 2D, if a simplicial mesh 
satisfies the following generalized Delaunay-type angle eondition (which is much weaker than that 
of equation (4.15 )): 


0 < 


fi fs-i + fi 




1 


+ -arccot 


+ -arccot 


det[iDf 

\ det[.Df 


cot(/3 ~-i) - 

^ pq,D • ' 


2e 


Q 5 


det[i:)Q^ 




2(t 




det[T)f^/] 


< TT 


(4.17) 


for every internal edge Cpq connecting the p-th and q-th vertices of the elements fig o,nd fig that 
share this common edge (see Figure 7), then we have the following three results: 

• The global stiffness matrix Kff is (reducibly/irreducibly) weakly diagonally dominant if 
a(x) > 0 and is (redueibly/irredueibly) strictly diagonally dominant if a{x) > 0 for all 
X G n. 
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The discrete single-field Galerkin formulation given by equations (2.24a)-(2.24b) in asso¬ 
ciation with equations (4.7)-(4.14) satisfies DwyVP k/TWMPk■ The quantities hp^cie 
are the heights of h and h are the heights of fi ~-i and /3_ 




are the relevant metric based dihedral angles in the elements Og and that face the edge 


Cpq, and the parameters 


I I loo,He 


“llooHe evaluated in Og based on the equa¬ 


tions (4.16a)-(4.16b). Similarly, ||t)|| _/ and ||a|| _/ are evaluated in TL . The constant 

oo,r2„ oo,r2p 


Qe equation (4.17) is given as follows: 


^q,Q.,n: ■= meas(Qg) 




00,f!e 


\a\ 


3h, 


+ 


00,f!e 


q,^e 


12 


+ meas(12g 


\v\ 


oo.Q^ 


a 


3h 


+ 






12 


(4.18) 


• Additionally, if the triangulation Th is interiorly connected, then the global stiffness matrix 
Kjf is irreducibly weakly or strictly diagonally dominant based on the nature of a{x) and 
the discrete single-field Galerkin formulation given by equations (2.24a)-(2.24b) satisfies 
DsMPk/DSMPk. 


Proof. For proof, see References [29,30]. 


□ 


Each method (global stiffness restriction method or local stiffness restriction method) has it 
own advantages and disadvantages. The advantage of the global stiffness restriction method is that 
we can operate at a global level. This gives us different types of relationships between various mesh 
parameters, which can be used in generating different types of triangulations, such as Delaunay- 
Voronoi, non-obtuse, well-centered, and anisotropic A4-uniform finite element meshes. For instance, 
Taylor series expansion of equations (4.15) and (4.17) gives the following restrictions on the metric 
based dihedral angles for all simplicial elements in a triangulation: 


^ ^ — 9 ^ (^ll'*^l|oo,Th + ^ Halloo, 77) 


(4.19a) 


0 < 




pq,D 




1 

-|—arccot 


/ 


det[T)j^ 
\ det[.Df 


cot(/3 




\ 

) 

/ 


1 


-|—arccot , . ~ 

2 I \ det|D„;. 


det[i:)!9, 


T I < ^ ^ {Hv\\oo,Th + /i^l|a||oo,rJ 


(4.19b) 


where the maximum element size h, maximum element normed velocity ||i?||oo, 77 ; and maximum 
element normed linear reaction coefficient ||a||oo ,77 are given as follows: 


h 

|'*^||oo,77 ■" 
l«lloo,r^ := 


■ max 

■ max 


a 


(4.20a) 

(4.20b) 

(4.20c) 


where /imax.fie is the maximum possible height in a given simplicial element fig, and 0{») is the 
standard “big-oh” notation [79]. Specifically, on a /i-refined triangular mesh, which conforms to 
Theorems 4.3 and 4.4, then equation (4.19a) implies that all the dihedral angles when measured in 
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the metric of have to be O {h\\v\\oo,Th + h‘^\\oi\\co,Th) acute/non-obtuse and equation (4.19b) 
indicates that the triangulation needs to be O {h\\v\\oo,Th + h‘^\\<^\\oo,Th) Delaunay. 

One downside of the global stiffness restriction approach is that obtaining mesh conditions is 
mathematically cumbersome. Moreover, extending it to non-simplicial low-order finite elements is 
extremely hard and is not straightforward. This is because the basis functions (pp^ spanning the 
hnite dimensional subsets and are multi-linear, which makes grad[(/?pg] on any arbitrary 
element fie to be non-constant. So most of properties given by equations (4.7)-(4.14) are not 
valid for low-order finite elements such as the Q4 element and its corresponding elements in higher 
dimensions. 

Now we shall describe the local stiffness restriction method and highlight the pros and cons of 
using this methodology. For the sake of illustration, we shall consider a pure anisotropic diffusion 
equation and assume that Xp = (0,0), Xg = (1,0), and Xr = (a, 6). Our objective is to find the 
coordinates {a,b), such that the local stiffness matrix is weakly diagonally dominant for any given 
type of diffusivity tensor. The local stiffness matrix for an anisotropic diffusion equation based on 
discrete single-field Galerkin formulation is given by: 


Ke= j BD{x)B^ dO (4.21) 

ne 

where the matrix B in terms of the coordinates (a, b) is given as follows: 

/ -6 (a - 1) \ 

B = -\ b -a (4.22) 

V 0 1 / 

Since the entries in the matrix B are constants, we have the local stiffness matrix given as 
follows: 


K, = B 


D{x) dO 


B^ 


(4.23) 


In the subsequent subsections, we present various sufficient conditions through which we can find 
these coordinates (a, h) and glean information on the possible angles and corresponding shape and 
size of the triangle PQR. 

4.2.2. T3 element for heterogeneous isotropic diffusivity. In this subsection, we consider 
the case where the diffusivity is isotropic and heterogeneous in the total domain. For this case, we 
show that the diffusivity does not have any influence on determining the coordinates {a,b). This 
means that the restrictions we obtain on the coordinates and the angles of the triangle PQR is 
independent of how the diffusivity is varying across the domain. The following is the local stiffness 
matrix for scalar heterogeneous isotropic diffusion: 

„ / b"^ + {a — 1)“^ a — a? — b"^ (a — I) \ 

Kf, = Jo D(a;) dll j a — — b^ a?' + b‘^ —a | (4.24) 

i \ (a-1) -a 1 / 
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where the integral of the diffusivity D{x) over the actual T3 element Og (triangle PQR) is given 
as follows: 


D = 


1 


meas(ne 


D{x) dn 


(4.25) 


where meas(ne) = | is the area of the T3 element PQR, which is always positive. The simplified 
expression for the local stiffness matrix is given as follows: 


K, = 


D 


6 " + (a-l)^ a-a^-¥ (a - 1 ) 


- -b"^ 


+ 6^ 
—a 


—a 

1 


(4.26) 


(a - 1 ) 

We shall now present the sufficient conditions so that the matrix Ke is weakly diagonally dominant: 


Condition No-1. Positive diagonal entries: {Kg)-- >0 Vi = 1,2,3. This restriction gives us 
the following inequalities: 


iKe)n = i - 1)') > 0 

{Kf,)22 +5)>0 

(-^e)33 ^ ^ ® 


(4.27a) 

(4.27b) 

(4.27c) 


As H > 0 and 6 > 0, it is evident that all of the inequalities given by equations (4.27a)-(4.27c) 
are trivially satished. Hence, this condition has no effect on obtaining restrictions on coordinates 
(a,Q. 


Condition No-2. Weak diagonal dominance of rows: | (iTe)ji | > 

i = 1,2,3 and j = 1,2,3. This restriction gives the following inequalities: 

(6^ + (a — 1)^) > (a^ + 6^ — a) + (1 — a) 

(yO? + 6^) > a + {a? + 6^ — a) 

1 > (1 — a) + a 


Vi,j, where 


(4.28a) 

(4.28b) 

(4.28c) 


Note that these inequalities (4.28a)-(4.28c) are trivially satished. Hence, this condition has no 
influence on obtaining restrictions on triangle PQR. 


Condition No-3. Non-positive off-diagonal entries: {Ke)^j <0 Vi 7 ^ j, where i = 1,2,3 and 
j = 1,2,3. This restriction gives the following inequalities: 




\Ke\2 

< (^e)i3 
^ (-^e)23 


= I (a - a 2 - 62 ) < 0 
= I (a - 1) < 0 Vi / i 

= i (-a) < 0 
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Using the fact that D > 0 and b > 0, and rearranging the above relations, we get the following 
inequalities: 



a < 1 
a > 0 


(4.29a) 

(4.29b) 

(4.29c) 


The region in which the coordinates (a, b) satisfy the above inequalities given by the equations 
(4.29a)-(4.29c) is shown in Figure 8. According to these inequalities (4.29a)-(4.29c), heterogeneity 
of the scalar diffusivity has no role in obtaining the feasible region for the coordinates {a,b). It is 
evident from Figure 8 that the interior angles of the triangle PQR are either acute or at most right- 
angle. Based on the sufficient conditions, one can also notice that an obtuse-angled triangle is not 
possible. So in order to satisfy discrete comparison principles, discrete maximum principles, and 
non-negative constraint, the triangulation of a given computational domain must contain acute- 
angled triangles or right-angled triangles. These three sufficient conditions show that non-obtuse or 
well-centered triangulations inherit all the three discrete versions of continuous properties of scalar 
heterogeneous isotropic diffusion equations. 

4.2.3. T3 element for heterogeneous anisotropic diffusivity. In this subsection, we con¬ 
sider the case where the diffusivity D(x) = ( ) is anisotropic and heterogeneous across 

\ Uxy\^) Uyy^'X.) J 

the domain. For the sake of brevity and ease of manipulations, we shall drop the symbol (x) in 
the components of the diffusivity tensor. Note that the symbol ‘x’ in the components of D(x) is 
dropped for the sake of convenience and should not be interpreted as though the diffusivity tensor 
is constant. As discussed in Section 2, the diffusivity tensor needs to satisfy certain properties. 
Based on equations (2.2) and (2.6), we derive various results related to D(x) that will be used in 
deriving mesh restrictions. 


Remark 4.5. In 2D, it is trivial to show that, if the matrix D{x) is symmetric, uniformly 
elliptic, and bounded above, then its components satisfy the following relations: 


Dxx ^ 0 

(4.30a) 

Dyy > 0 

(4.30b) 

Dxx Dyy > D,j.y 

(4.30c) 


Let us denote e := -2^ and y := -2^, where D^x, Dxy, and Dyy are the components of matrix 

j-^xx J-^XX 

given by equation (4.7). From Theorem 4.2, it is evident that Dxx > 0, Dyy > 0, and 
DxxDyy > ^xy So from equation (4.30c), we have rj G These two non-dimensional 

quantities e and y govern the mesh restrictions that we impose on the coordinates {a,b). From 
equations (4.21) and (4.22), the stiffness matrix for any given anisotropic diffusivity tensor is given 
as follows: 


Ke 


„ 2 b „ 

D^P3^^-D^y(h—’lah)^-Dyya{a—\) 

2b 

Dxy^~\~Dyy{_^ 1 ) 

\ 2b 


Dxxb ~\-Dxy(^b 2(lb'j-\-DyyCL(^CL l) 

I I 

Dxxb 2Dxy(lb-\- Dyyd 

I Wi 

Dxyb Dyyd 
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Dxyb-\-Dyy{d 1 ) 

I 261 

Dxyb Dyyd 

3^ 

^yy 

26 


\ 

/ 


(4.31) 


We now present sufficient conditions so that the matrix K^, is weakly diagonally dominant. 
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Condition No-4. Positive diagonal entries: {Kg)-- >0 Vi = 1,2,3, gives the following 
relations: 


iKe)u = 


iKe)n 

(K, 


_ Dxxb^—2Dxyb(a—l)+Dyy(a—l)^ 


2 b_ 

Dxxb ‘2DxyOib-\-DyyOi 


> 0 


■e )22 


2b 


> 0 


(K. 


eJ33 


_ ^yy ^ n 

~ 2b ^ 


As Dyy > 0 and 6 > 0, rearranging the above relations, we have the following restrictions: 

{Ke)ii = (^b\/-\a- l\^JDyy'^ + 26|a - 1| Dy,y,\Jbyy - Sgii [{tt - l|] ^ > 0 (4.32a) 

{Ke )22 = {b\l Dxx - \a\\lby^ 2h\a\ (^\f^x\[byy - Sgn [\a\] > 0 (4.32b) 

(Ke )33 = ^ > 0 (4.32c) 

where Sgn[«] is the standard lignum function (which provides the sign of the real number). It is 

evident that ^D^x\JDyy > Dxy Hence, equations (4.32a)-(4.32c) are trivially satisfied for any 
abscissa a. 


Condition No-5. Non-positive off-diagonal entries: (Kg)- <0 Vi 7 ^ j, where i = 1,2,3, 
and j = 1,2,3. This restriction gives the following relations: 




^ (^e)i3 
, {^e)23 


_ DxX^ ‘2Q,b'j-\-]Dyy(l(^(l 1 ) Q 

^ ^ 2b — 

_ —Dxyb^^yy{<^—^) ^ 

“ ^ Jb 

_ Dxyb—Dyya ^ p. 

~ 2b — ^ 


Using the parameters e, 77 , and the fact that ordinate 6 > 0, we have the following inequalities: 



Q - 1 ^ V 
b - e 

a > V 

b - e 


(4.33a) 

(4.33b) 

(4.33c) 


which dictate the feasible region for coordinates (a, 6 ). For a given e and by varying rj, which lies 
between —^/e and ^/e, we get different feasible regions for (a, 6 ). Herein, we have chosen e = 10 
and r] S {—1,0,1}. For these values, we have plotted the feasible region based on the inequalities 
(4.33a)-(4.33c). From Figures 9-14, the following can be inferred based on the feasible region: 

• If T/ = 0, the possible T3 elements are either acute-angled or right-angled triangles. 

• If either r] < 0 or rj > 0, then obtuse-angled triangles are also possible. Moreover, the 
resulting triangles can be skinny or skewed. 
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Condition No-6. Weak diagonal dominance of rows: | (-FCe)ji 
i = 1, 2,3, and j = 1, 2, 3. This gives the following relations: 


> EI 


yi,j, where 


(6^ — 2776(0 — 1) + e(a — 1 )^) > (6^ + 77(6 — 2ab) + ea{a — 1)) + (776 — e(a — 1 )) (4.34a) 

(6^ — 27706 + eo^) > (6^ + 77(6 — 2 a 6 ) + eo(a — 1)) (eo — 776) (4.34b) 

e > (776 — e(a — 1)) + (eo — 776) (4.34c) 


if Condition-4 and Condition-5 are satisfied, then this condition is trivially satisfied. In a similar 
fashion, for a general case, where Xp = (xi,7/i), Xg = {x 2 ,y 2 ), and Xr = ( 2 : 3 , 7 / 3 ), we have the 
following conditions based on the local stiffness restriction method: 


(yi - 7/3)(ys - y2) - v{xi - Xs){y 3 - ^2) - 77(3:3 - X 2 ){yi - ys) + e(xi - 0:3)(X3 - 2:2) < 0 ( 4 . 35 a) 

(y2 - yi)(y3 - ^2) - 77(^3 - X2)(y2 - yi) - 77(3:2 - Xi){y3 - ^2) + e{x2 - 3:i)(x3 - X2) < 0 ( 4 . 35 b) 

(yi - y3)(y2 - yi) - 77(^1 - X3)(y2 - yi) - 77(3:2 - 3:i)(yi - ys) e(xi - 3:3 )(x 2 - 3:1) < 0 ( 4 . 35 c) 

( 3:1 - X3)(7/2 - 7 / 1 ) - {xi - X2)(y3 - yi) > 0 (4.35d) 


where the hrst three inequalities given by equations (4.35a)-(4.35c) are obtained based on the 
condition that {Ke)-j < 0. The last inequality given by equation (4.35d) is the result of the 
condition that meas(fle) > 0- 

The beneht (attractive feature) of the local stiffness restriction method is that the local stiffness 
matrix for the discrete Galerkin formulation given by equations (2.24a)-(2.24b) can be calculated 
quite easily and could be extended to even non-simplicial elements (see subsection 4.4). Using this 
approach, we can obtain general restrictions and analytical expressions relating various coordinates 
of an arbitrary mesh element, using popular symbolic packages like Mathematica [80]. But a flip-side 
of this procedure is that incorporating the inequalities given by equations (4.35a)-(4.35d) in a mesh 
generator is very difficult and needs further detailed investigation. Additionally, the conditions 
obtained using this method are stringent and similar to that of Theorem 4.3 given by the global 
stiffness restriction method (which will be evident based on the numerical examples discussed in 
the following subsection). Finally, it should be noted that extending the local stiffness restriction 
method to include advection and linear reaction is straightforward and shall not be dealt with to 
save space. We shall now present various numerical examples and respective triangular meshes 
corresponding to different types of D(x). Using these meshes, we shall analyze and study in detail, 
which kind of DMPs and DCPs are preserved. 


4.3. Numerical examples based on different types of triangulations. In this subsec¬ 
tion, we shall first briefly discuss on a metric tensor Ai{x) to satisfy DCPs, DMPs, and NC. Based 
on this metric tensor, we shall describe an algorithm to generate various types of DMP-based tri¬ 
angulations (mainly utilizing open source mesh generators, such as Gmsh and BAMG). Simplicial 
meshes constructed based on Ai{x) (where the metric tensor Ai{x) is not equal to a scalar multiple 
of identity tensor) are called anisotropic A4-uniform simplicial meshes. They are uniform in the 
metric specified by Ai{x) [61,72] and are of primal importance in satisfying various important 
discrete properties in the areas of transport of chemical species, fluid mechanics, and porous me¬ 
dia applications [81-83]. Given a Ai{x), there are different approaches to generate anisotropic 
A4-uniform simplicial meshes. Some of the notable research works in this direction include blue 
refinement [84], bubble packing [85], Delaunay-type triangulation [81], directional refinement [86], 
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front advancing [85], local refinement and modification [87], and variational mesh generation [88]. 
In general, the metric tensor Ai{x) is symmetric and positive definite, and gives relevant informa¬ 
tion on the shape, size, and orientation of mesh elements in the computational domain. 

Let X be an affine mapping from the reference element flref to background mesh element fig. 
Denote its Jacobian by Jn^. The affine mapping X and its Jacobian are given as follows: 

a: = X(^, N) = (4.36a) 

= X^DN (4.36b) 

where X is the nodal matrix, which comprises of nodal vertices {xi, X 2 , ■ ■ ■, a^nd+i) of an arbitrary 
simplex Dg. The elements of vector N consists of shape functions corresponding to Di-ef. The 
entries of matrix DN (which are constants for simplicial elements) correspond to the derivatives of 
shape functions. Straightforward manipulations on equations (4.8) and (4.36b) give the following 
result: 


9p,ref (4.37) 

where Qp j-ef is the corresponding q-vector of the reference element Di-ef. Huang [72] has shown that 
an anisotropic A4-uniform simplicial mesh satisfies the following two conditions: 

p^^meas{Qe) = VHg G % (4.38a) 

^tr [JlMn^Jn.] = det [JlMnJn.] VHg G Th (4.38b) 

The quantities given as follows: 

Mn^ :=-^ [M{x)dn VHg G Th (4.39a) 

meas(i2e) J 

He 

Pn^ '■= \/det [Mn^ VDg G % (4.39b) 

(Jh ■■= E p^jneas{^e) (4.39c) 

fle&Th 

The condition given by the equation (4.38a) is called equidistribution condition and that by equation 
(4.38b) is called alignment condition. Equidistribution condition decides the size of Dg, while 
alignment condition characterizes the shape and orientation of Dg. From AM-GM inequality [89], 
equation (4.38b) implies the following: 

= I VHgGTft ( 4 . 40 ) 

Now, through trivial manipulations on equations (4.37), (4.10), and (4.40), the metric tensor 
has to satisfy the following equation in order to meet various DMPs: 

Mn. = QnAl (4-41) 

where is an arbitrary piecewise positive scalar constant, which in general is a user-define 
parameters. Loosely speaking, an anisotropic A4-uniform simplicial mesh satisfying (4.41) satisfies 
weak DMPs. Furthermore, if the resulting simplicial mesh is interiorly connected, then it satisfies 
strong DMPs. 
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For a given type of metric tensor (for example, Ai{x) given by equation (4.41)), open source 
mesh generators (such as BAMG, BL2D, and MmgSd) take this as an input and operate on a 
background mesh to produce an anisotropic A4-uniform simplicial mesh. Algorithm 1 provides a 
methodology to develop an anisotropic A4-uniform simplicial mesh based on a background mesh, 
and Figure 10 highlights the salient aspects of this algorithm. Nevertheless, it should be noted 
that Algorithm 1 is a general algorithm to generate DMP-based triangulations (is not limited to 
either Gmsh or BAMG) for any type of mesh generator, which operates on a background mesh. On 
the other hand, there are certain open source and commercially available software packages, such 
as CGAL [90] and Simmetrix [91], which create an anisotropic A4-uniform simplicial mesh directly 
based on the metric tensor [92,93] without the need of background meshes. Investigation of such 
mesh generators is beyond the scope of this research paper and is neither critical nor central to 
the ideas discussed here. Before we discuss various numerical examples based on different types 
of triangulations, we shall now present certain important (mesh-based) non-dimensional numbers 
relevant to our numerical study. Such a discussion is first of its kind and is not discussed elsewhere. 

Remark 4.6. It should he noted that STEP-7 in Algorithm 1 might he computationally intensive, 
as we need to solve a series of small constrained optimization problem for each £ Th,i ’ and at 
every iteration level ‘i The corresponding constrained optimization problems are given as follows: 

VI < z < Maxiters (4.42a) 

VI < i < Maxiters (4.42b) 


u||e,i ;= maximize ||u(a;)|| 
a||e,i := maximize a{x) 


As is convex, triangular, and closed, by half-space representation theorem for convex polytopes 
[95, Section-2.2.4] equation (4.42a) can he written as follows: 

maximize ||u(a3)|| (4.43a) 

subject to Ae^iX ^ VI < i < Maxiters (4.43b) 

where is a 3 x 2 matrix and is a 3 x 1 vector, whose coefficients correspond to the linear 
inequalities defining the relevant half-spaces and supporting hyper-planes of the triangle If the 
element maximum value ‘||u||e,i ’ is known a priori (through analytically or by means of a rigorous 
mathematical analysis), then one can use such information in STEP-7 of Algorithm 1. Otherwise, 
we need to solve equations (4.43a)-(4.43b) using the standard constrained optimization algorithms 
for small-scale problems [95]. Similarly, equation (4.42b) can he reformulated based on the lines of 
equations (4.43a)-(4.43b). 

4.3.1. Peclet and Damkohler numbers for simplicial meshes. Herein, we shall describe 
three types of Peclet and Damkohler numbers for simplicial meshes. They are devised based on 
Theorems 4.3 and 4.4. However, it should be noted that extending it to non-simplicial elements, 
such as Q4, is not straightforward. This is because in order to construct Peclet and Damkohler 
numbers for non-simplicial meshes, one needs to obtain mesh restrictions using the global stiffness 
restriction method. This is beyond the scope of the current paper. 
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Algorithm 1 An iterative method to generate an anisotropic Al-uniform mesh satisfying discrete 
principles 

1: INPUT: Background mesh (7h,0) NeleQ, Nvq, and Nbvo)', anisotropic diffusivity tensor (D(x)); 
velocity vector field (v(x)); linear reaction coefficient (a(x)); maximum number of iterations 
(Maxiters); piecewise positive scalar element metric constants a stopping 

criteria (StopCrit) 

• Th,o is the initial background triangulation on which an anisotropic mesh generator 
operates 

• Nvq and NbvQ are correspondingly the total number of vertices and boundary vertices 


2 

3 

4 


5 : 

6 : 


7 : 


8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 


Set the iteration number: i = 0 
while (True) do 

Compute the element average anisotropic diffusivity tensor using a quadrature rule (for ex¬ 
ample, see References [94]) 


i := 


D(x)dQ Ve = 1,2, ■■■ jNelci 


' meas(fle) 

Compute the element metric tensor by explicitly inverting Df,,i 

• Me,i ■■= &e{De,i)~^ Me = 1,2, ■■■, Nelci 

Based on the set of metric tensors compute a new triangulation Th,i 

• Output the new triangulation Th^i- Corresponding to this Th^i, we have Nelei, Nvi, 
and Nbvi 

Compute the following quantities: Ve = 1, 2, • • • ,Nelei 


Dei] K . 

’ miTi 


\V\ 




I 00,Qe i 


• Need to use a constrained optimization methodology to calculate ||ii||e,i and ||a||e,i 
(see Remark 4.6 for more details) 

if (StopCrit = Anisotropic non-obtuse angle condition)then 

Check the inequality given by equation (4.15) in Theorem 4.3 Ve = 1, 2, • • • , Nelei 
if (true) then 

OUTPUT: The triangulation Th,i and corresponding {M.e,i}e=iM EXIT 
else 

Update Th,i i— Th,i, Nelei ■(— Nelei, Nvi ■(— Nvi, Nbvi •(— Nbvi, and z •(— (z + 1) 

end if 
end if 

if (StopCrit = Generalized Delaunay-type angle condition)then 

Check the inequality given by equation (4.17) in Theorem 4.4 Ve = 1,2, • • • , Nelei 
if (true) then 

OUTPUT: The triangulation Th,i and corresponding {M.e,i}e=iM EXIT 
else 

Update Th^i i— Th^i, Nelei •(— Nelei, Nvi •(— Nvi, Nbvi ■(— Nbvi, and z •(— (z + 1) 

end if 
end if 

if (z > Maxiters) then 

OUTPUT: The existing triangulation Th,i and corresponding 
Anisotropic A4-uniform triangulation not found in Maxiters. EXIT 
end if 
end while 








• Element Peclet and Damkohler numbers: Based on Theorem 4.3, one can define the follow¬ 
ing mesh-based non-dimensional element Peclet (Peo^) and Damkohler (Dan^) numbers: 


Peoe 




^max.He ll'*^llcxD,r2e 

min ^ DQ ^ 

^maXjOe ^pumaXjOe ll*a||gQ 


A 


min,Dn, 


where the height /ipumax.fie given as follows: 


(4.44a) 

(4.44b) 


0 < hi < /i 2 “El ''' “El hi “El ''' “El ^pumax,f2e — ^max,r 2 e — f; 2, • • • , nd -\- 1, Dg £ 7/i (4.45) 


Correspondingly, using equations (4.44a)-(4.44b) in equation (4.15) gives the following 
(stronger) mesh restriction condition based on Theorem 4.3: 


0 < 


Peoe 


{nd+l) cos(/3 . ~-i; 

J > Qa 


+ 




(nd+l) {nd + 2) cos(/3 .--y 

J ! r 2 e 


< 1 


i = max and j = pumax, i ^ j, VDg G Th 


(4.46) 


• Edge Peclet and Damkohler numbers: In a similar fashion, utilizing Theorem 4.4, one can 
define the following mesh-based non-dimensional edge Peclet (Peo^^ep,) and Damkohler 
(DaQ^,epJ numbers: 


Pe^e 


■,epq 




meas(De)||^||^^^ 

meas(De)||a|loo,T7, 

\Jdei[bnJ 


(4.47a) 

(4.47b) 


Correspondingly, using equations (4.47a)-(4.47b) in equation (4.17) gives the following 
(weaker) mesh restriction condition based on Theorem 4.4: 


0 < 


27r 




( 


pqDc 


H-arccot 

27r 


1 


V 


\ det[bn,' 


cot(B ~-i) — 
pq,D^' 


2Pe, 


q' f. 


^q' p 




6 


H-arccot , . ~ 

21 I \ det|D„, 


ho,(,3 _ j_2P^ 

^^pq,Dn, ^ 3 




2Pe, 


Og,epg 


Ba, 


Q^,epq 


< 1 


(4.48) 


• Global mesh Peclet and Damkohler numbers: On sufficiently fine /i-refined simplicial 
meshes (which confirm to Theorem 4.3), one can define a global mesh Peclet (Pe/i) and 
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Damkohler (Da/i) numbers by modifying equations (4.44a)-(4.44b) as follows: 


Feh := 

Ba/j := 


h max 
ne&Th 



min 


A 


min,DQ 


e 


/i^ max 
^eG'Th 



min 

^e^Th 


A 


min,DQ 


e 


(4.49a) 


(4.49b) 


Conservatively, equation (4.46) can be modified to give a (stronger) global mesh restriction 
condition based on Theorem 4.3: 


0 < 




(nd + 1) min 
HeeTh 


T + 


Ba/i 


cos{B 


{nd + 1) {nd + 2) min 


cos{B ~-i) 


< 1 


i = max and j = pumax, i ^ j 


(4.50) 


One should note that equation (4.50) can provide a useful a priori (conservative) estimate 
on h for constructing highly refined simplicial meshes. 

4.3.2. Physics-based Peclet and Damkohler numbers. For isotropic diffusivity, it is well- 
known that the following three physics-based non-dimensional numbers can be used to understand 
the qualitative nature of the solutions [5,96,97]: 

max [||v(x)||oo] L 

Pe£) := ,, - Peclet number (4.51a) 

mm [i4(x)J 

max[||a(x)||oo] L 

Bai := —j-—-—7—— Damkohler number of first kind (4.51b) 

max[||v(x)||ooJ 

tcgf! 

max[||a(x)||oo] 

Ban D ■— -r—- Damkohler number of second kind (4.51c) 

’ min [D(x)J 

where L is the characteristic length of the domain and || • jjoo is the standard max-norm/infinity- 
norm for vectors [56]. However, it should be noted that the above three non-dimensional numbers 
are not independent, as they satisfy the relation Ban = Pez) Baj. Physically, the non-dimensional 
Peclet number characterizes the relative dominance of advection as compared to diffusion processes. 
For larger Peclet numbers, the advection process dominates and for smaller Peclet numbers, the 
diffusion process dominates. 

In the literature on chemically reacting systems (e.g., see [98, Table 22.2.1]), there exists various 
Damkohler numbers that relate progress of chemical reactions with respect to mixing, diffusivity, 
reaction coefficient, thermal effects, and advection. The Damkohler number of first-kind, Bai, gives 
information about the relative influence of a linear reaction coefficient to that of advection. For 
small values of Baj, advection progresses much faster than decay of the chemical species and has 
the opposite effect for large values of the number. The Damkohler number of second kind. Ban 
gives information related to the progress of chemical reaction with respect to diffusion. Based on 
the chemical system under consideration, these three non-dimensional numbers dictate how the 

33 





















reaction, advection, and diffusion interact with each other. On the other hand, one should note 
that extending it to anisotropic diffusivity is not straightforward. Motivated by equations (4.49a)- 
(4.49b), a way to define physics-based non-dimensional numbers for anisotropic diffusivity tensor 
D (x) is as follows: 

max[||v(x)||oo].b 

PeD := r - y Peclet number (4.52a) 

min J 

xGQ 

max [||a(x)||oo] 

D —FT-Damkohler number of second kind (4.52b) 

min J 


where is the minimum eigenvalue of D(x) at a given point x. The Damkohler number 

of first kind, Dai, given by equation (4.51b) remains the same as it does not depend on diffusivity 
tensor D(x). Alternatively, inspired by equations (4.47a)-(4.47b), one can define a different set of 
physics-based Peclet and Damkohler numbers. These are given as follows; 


PeD ;= 


Daii^D : 


max[||v(x)||oo] L 

xGQ 



max j 

x£Q 


max[||a(x)||oo]T^ 

xGQ 



m^ [A 'max,D(x)] 


Peclet number (4.53a) 


Damkohler number of second kind (4.53b) 


where ^rnax,T>{-x.) is fh® maximum eigenvalue of D(x) at a given point x. One should note that both 
these sets of non-dimensional numbers (given by equations (4.52a)-(4.52b) and (4.53a)-(4.53b)) 
are perfectly valid, as anisotropy in diffusivity tensor can introduce multiple ways of defining 
physics-based Peclet and Damkohler numbers. Now, we shall present various numerical examples 
to demonstrate the pros and cons of the mesh restrictions approach. 

4.3.3. Test problem #1: Transport in fractured media. This test problem has profound 
impact in simulating the transport of chemical species in fractured media [99]. The numerical 
simulations, using the Delaunay-Voronoi triangulation (with Maxiters = 50) for different cases of 
diffusivity, velocity field, and linear reaction coefficient, are presented in Figure 15. Homogeneous 
Dirichlet boundary conditions are prescribed on the sides of the fractured domain; cP(x) = 1.0 
on the left set of fracture lines and cP(x) = 2.5 on the right set of fracture lines. The volumetric 
source /(x) is zero inside the fractured domain. Diffusivity is assumed to isotropic and scalar, 
whose value is given by D(x) = 10“^. We perform numerical simulations for three different cases 
of velocity field and linear reaction coefficient, which are given by v(x) = (0.0,0.0) and a(x) = 0.0, 
v(x) = (0.1,1.0) and q;(x) = 0.0, and v(x) = (0.1,1.0) and a(x) = 1.0. 

From Figure 15, it can be inferred that we need highly refined DMP-based triangular meshes to 
obtain physically meaningful values for concentration for large values of edge Peclet and Damkohler 
numbers. The white region in the figures indicates the area in which the value of concentration is 
negative and also violated the maximum constraint. The coarse Delaunay-Voronoi mesh obtained 
using the open source mesh generator Gmsh satisfies NC and DMPs in the case of pure diffusion. 
But, this is not true for AD and ADR cases. In such scenarios, it produces unphysical values for 
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the concentration field. Moreover, the percentage of nodes that have violated NC and maximum 
constraint is also very high. 

Quantitatively, Tables 1 and 2 provide more details pertinent to the violations in NC and DMPs 
for AD and ADR cases. It should be noted that the decrease in unphysical values of concentration 
is not monotonic. This is because the generalized Delaunay-type angle condition, in general, 
does not ensure uniform convergence for diffusion-type equations in L°° norm (for example, in case 
of pure isotropic diffusion, see the mesh restriction result by Ciarlet and Raviart [12]). Figure 16 
shows that the weak DMP-based condition is satisfied only for pure diffusion equation. But in all 
other cases, generalized Delaunay-type condition is violated. In the case of pure diffusion, it 
should be noted that DMP-based mesh given in Figure 16 is not interiorly connected. Hence, it 
only satisfies DWMP^, but not DSMP^- 

4.3.4. Test problem #2: Species dispersion in subsurface flows. A pictorial description 
of the boundary value problem with various parameters is shown in Figure 17. Relevant (coarse) 
background mesh used and the corresponding DMP-based mesh (with Maxiters = 50) obtained 
using BAMG are shown in Figure 18. The diffusivity tensor for this problem is taken from the 
subsurface hydrology literature [ 100 ] and is given as follows: 

Dsubsurface(x) = aT||v||H-- V (g) V (4.54) 

where (g) is the tensor product, I is the identity tensor, v is velocity vector field of the subsurface flow, 
and ot and are, respectively, transverse and longitudinal diffusivity coefficients with ot = 0.01 
and aL = 0.1. It should be emphasized that we have neglected advection. Correspondingly, 
the numerical values for the velocity vector field used to define the diffusion tensor is given by 
v(x) = (1.0,1.0). This test problem has importance in simulating diffusion of chemical species in 
subsurface flows of hydrogeological systems [101-103]. 

Numerical simulations using these meshes are shown in Figure 19. The white region in this 
figure depicts the area in which the value of concentration is negative. From Figure 19, it is 
apparent that the coarse anisotropic triangulation violates the DMPs and NC. This is because the 
Algorithm 1 did not converge in Maxiters. However, quantitatively, this violation in NC is low 
as compared to background mesh. Specifically, the minimum concentration and the percentage of 
nodes that have violated the non-negative constraint on the background mesh is about —4.8 x 10“^ 
and 2.13%, while these values on the anisotropic triangulation are around —1.35 x 10~® and 0.34%. 
Additionally, from Figure 19, it is evident that we need a highly refined DMP-based anisotropic 
mesh to avoid negative values for concentration. Nevertheless, it should be noted that there is a 
considerable decrease in the negative values for concentration if a traditionally h-rehned anisotropic 
triangulation is used (see Figure 23 and subsubsection 4.3.6 for more details). 

4.3.5. Test problem #3: Contaminant transport in leaky wells. The computational 
domain is a circle with a hole centered at origin (0,0). The radius of the circular hole and the 
circular domain are 0.1 and 1.0. Numerical simulations are performed for four different cases of 
the velocity vector field and linear reaction coefficient, which are given by v(x) = (0.0,0.0) and 
a(x) = 0.0, v(x) = (1.5,1.0) and a(x) = 1.0, v(x) = (5.0,0.5) and q;(x) = 1.0, and v(x) = (0.0,0.0) 
and q:(x) = 1000. Each case is designed to test a particular aspect. For example, v(x) = (5.0,0.5) 
and q;(x) = 1.0 corresponds to advection-dominated ADR problems, while v(x) = (0.0,0.0) and 
a(x) = 1000 corresponds to reaction-dominated diffusion-reaction problems. The diffusivity tensor 
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for this problem is given by equations (1.1)-(1.3). Correspondingly, the values for the parameters 
c^rnax, c^min, and 6 are equal to 1, 0.001, and 

The computational domain in this test problem has various practical applications related to well 
design in multi-aquifers and groundwater distribution systems [104 106], understanding emanation 
of gaseous hydrocarbons through bore-holes in oil and gas reservoirs [107 109], and to study 
leakage of contaminants (such as CO 2 , salts, and nitrates) through abandoned wells [110-113]. 
The background mesh used and the corresponding DMP-based (coarse) mesh obtained using BAMG 
are shown in Figure 20 (Maxiters = 50). For the cases when v(x) = (0.0,0.0) and a(x) = 0.0, 
v(x) = (1.5,1.0) and a(x) = 1.0, Algorithm 1 converged in Maxiters. But for v(x) = (5.0,0.5) 
and a(x) = 1.0, v(x) = (0.0,0.0) and a(x) = 1000, Algorithm 1 did not converge in Maxiters. 
Herein, the DMP-based coarse mesh is composed of needle-type triangles. This is because the ratio 
of the minimum eigenvalue of D(x) to its maximum is 0.001 (which is related to the aspect ratio 
of the sides of the triangle in the DMP-based anisotropic mesh). Moreover, it is evident that the 
triangles in the mesh are aligned and oriented along the principal axis of the eigenvectors of the 
diffusivity tensor. 

Numerical simulations using these meshes are shown in Figure 21. The white region in the 
figures (circular annulus) shows the area in which the value of concentration is negative. The mini¬ 
mum concentration and the percentage of nodes that have violated the non-negative constraint for 
the background mesh are about —1.67 x 10“^ and 30.28%. As the anisotropic mesh is coarse and the 
Algorithm 1 did not converge in Maxiters for advection-dominated ADR problems and reaction- 
dominated diffusion-reaction problems, the resulting mesh not only violates the non-negative con¬ 
straint, but also produces spurious oscillations. The minimum concentration and the percentage of 
nodes that have violated the non-negative constraint for the case when v = (5.0,0.5) and a = 1.0 
are about —1.78 x 10“^ and 13.47%, whereas for the case when v = (0.0, 0.0) and a = 1000 are 
around —2.79 x 10“^ and 20.54%. Hence in both these cases, we need a highly refined DMP-based 
mesh to avoid negative values for concentration (see subsubsection 4.3.6 for more details). For 
scenarios when v(x) = (0.0,0.0) and a(x) = 0.0, v(x) = (1.5,1.0) and a(x) = 1.0, the coarse 
DMP-based mesh is interiorly connected and satisfies DsMP^ (or DSMP/^, when a(x) = 0). 

4.3.6. Issues with DMP-based h-refinement. From the above set of test problems, it is 
apparent that mesh refinement is needed to avoid spurious node-to-node oscillations and satis¬ 
faction of various discrete principles (mainly for the cases when the values of the velocity vector 
field and linear reaction coefficient are predominant as compared to minimum eigenvalue of dif¬ 
fusivity tensor). In general, within the context of computational geometry and mesh generation 
literature, there are various methods to generate different types of /i-refined meshes [61,63,73]. 
However, it is evident from these test problems that we are interested in /i-refined meshes that 
conform to either anisotropic non-obtuse angle condition or generalized Delaunay-type 
angle condition. Herein, we shall study two different and important methodologies that are 
popular in mesh generation literature and implemented in the open source mesh generators, such 
as Gmsh and BAMG in FreeFem-|— |-. Our objective is to understand whether the h-refined meshes 
generated using these mesh generators satisfy DMPs, DCPs, and NC or not. 

Traditional h-refinement. In this method, an h-rehned mesh is obtained by splitting each tri¬ 
angle in the coarse mesh into multiple triangles. Based on this methodology, numerical simulations 
are performed using the traditional h-refined meshes (which are derived based on the DMP-based 
coarse meshes presented in Figures 15, 18, and 20). The resulting concentration profiles for test 
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Table 1. Fractured domain with isotropic diffusivity: For AD equation 


Delaunay triangulation 
{Nv,Nele,h) 

Concentration 

% of nodes violated 

Minimum 

Maximum 

(Non-negative constraint) 

(Maximum constraint) 

(539,906,0.1) 

-10.51 

9.85 

17.44 

13.73 

(1564,2826,0.05) 

-18.26 

5.86 

18.22 

14.19 

(5090,9620,0.025) 

-4.65 

5.71 

16.09 

13.77 

(18372,35665,0.0125) 

-2.97 

5.17 

12.82 

12.70 

(69995,137881,0.00625) 

-1.62 

3.89 

8.08 

8.47 


problems #1, ^2, and ^3, are shown in Figures 22-24. Qualitatively and quantitatively, the 
following inferences can be drawn from these numerical results: 

• Test problem ^1: For isotropic diffusivity, from Figure 22, Table 1, and Table 2, it is 
apparent that there is a decrease in negative values for concentration and reduction in 
spurious oscillations. 

• Test problem ^2: Qualitatively, based on Figure 23, it can be concluded that there is 
a considerable decrease in the violation of non-negative constraint. Quantitatively, for 
this /i-refined anisotropic triangulation, the minimum concentration and the percentage of 
nodes that have violated the non-negative constraint is about —1.15 x 10“^^ and 0.02% 
(which is significantly close to machine epsilon Cmach ~ 2.22 x 10“^®). 

• Test problem ^3: For the pure anisotropic diffusion case, the coarse anisotropic DMP- 
based mesh given in Figure 20 satisfies NC and all of the DMPs. However, contrary to this, 
on traditional h-refinement, the resulting anisotropic triangulation produces unphysical 
negative values and violates the DMPs. Quantitatively, this violation is far from machine 
epsilon, emach- Correspondingly, the minimum concentration and the percentage of nodes 
that have violated the non-negative constraint is about —2.7 x 10“^ and 28.51%, which is 
way higher as compared the numerical simulations based on the background mesh. 

To summarize, it is clear that traditional h-refinement does not always reduce the unphysical nu¬ 
merical values. Certainly, there is a decrease in negative values for concentration for test problem 
#1 and ^2. But, this is not the case for test problem ^3. This is because the methodology to 
obtain traditional h-refinement meshes from a given coarse mesh need not conform to the conditions 
outlined in either Theorem 4.3 or Theorem 4.4. 

Non-traditional h-refinement. In this method, a h-rehned mesh is obtained directly from 
the background mesh (using Algorithm 1) by change certain parameters related to metric tensor, 
geometry of the domain, and number of nodes on the boundary of the domain [ 20 , Chapter-5]. 
It should be noted that this methodology is completely different from the traditional h-rehnement 
procedure, as we never generate a coarse DMP-based triangulation and then split the respective 
mesh elements. However, even this non-traditional h-refined mesh is not guaranteed to produce 
physics-compatible numerical values for concentration (as the mesh generation procedure need not 
converge in Maxiters). This is evident from the numerical simulations performed on the non- 
traditional /i-refined mesh for test problem f^3. The corresponding concentration profile is shown 
in Figure 25. Quantitatively, the minimum concentration and the percentage of nodes that have 
violated the non-negative constraint are about —9.69 x 10“^ and 34.11%, which is comparatively 
similar to that of the traditional h-refinement methodology. 
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Table 2. Fractured domain with isotropic diffusivity: For ADR equation 


Delaunay triangulation 
{Nv,Nele,h) 

Concentration 

% of nodes violated 

Minimum 

Maximum 

(Non-negative constraint) 

(Maximum constraint) 

(539,906,0.1) 

-9.64 

4.53 

14.29 

8.16 

(1564,2826,0.05) 

-1.54 

3.57 

18.03 

1.53 

(5090,9620,0.025) 

-4.47 

3.30 

16.01 

0.29 

(18372,35665,0.0125) 

-2.95 

2.56 

12.89 

0.04 

(69995,137881,0.00625) 

-1.62 

2.50 

8.06 

0.00 


Table 3. Errors in local and global species balance: For pure isotropic and anisotropic 
diffusion equation. 


h 

Test problem jfl 

h 

Test problem (f2 

Local error (abs. max.) 

Global error 

Local error (abs. max.) 

Global error 

0.1 

1.80 X 10"3 

1.53 X 10"2 

0.08 

5.75 X 10-^ 

-1.44 X 10-^ 

0.05 

4.66 X 10“'’^ 

1.23 X 10"2 

0.042 

1.47 X 10-'’^ 

-1.38 X 10-^ 

0.025 

3.48 X lO""^ 

8.48 X 10"3 

0.028 

7.34 X 10-3 

-1.12 X 10-^ 

0.0125 

7.14 X 10“^ 

7.82 X 10-3 

0.0135 

1.59 X 10-3 

-6.83 X 10-3 

0.00625 

6.75 X 10"^ 

5.38 X 10-3 

0.0075 

5.46 X 10-3 

-4.36 X 10-3 


4.3.7. Errors in local and global species balance. It is well-known that without using 
a post-processing method, the discrete standard single-field Galerkin formulation does not possess 
local and global conservation properties [46,114]. In finite element literature, there are various post¬ 
processing methods to quantify the errors incurred in satisfying local and global species balance [5, 
96,115,116]. Herein, we shall use Herrmann’s method of optimal sampling [116, Subsection 15.2.2], 
which is a popular post-processing technique to obtain derived quantities from primary variables 
(such as the concentration variable in single-field finite element formulations). In the context 
of recovery-based error estimators [62], this post-processing method is also known as traditional 
global smoothing method [115, Subsection 9.9]. To summarize, this method involves minimizing 
a unconstrained least-squares flux functional to obtain nodal flux vectors. Then, using these flux 
vectors, local and global species balance errors are computed. 

Qualitatively, the contours corresponding to local species balance errors for test problems 
^2, and 7)^8 on coarse DMP-based meshes are shown in Figure 26. From this figure, it is clear that 
test problem ^3 exhibits considerable errors in satisfying local species balance. Quantitatively, 
Table 3 provides local and global species balance errors on traditional /i-refined meshes for test 
problems ^1 and ^2. From this table, it can be concluded that the decrease in these errors is slow. 
Hence, there is need for locally conservative DMP-preserving finite element methods. 

Remark 4.7. It should be noted that there are various ways in which one can develop lo¬ 
cally conservative global smoothing methods. For example, this can be achieved by constructing a 
constrained monotonic regression based method (following the procedure outlined by Burdakov et 
ah [46]) to obtain a conservative flux vector. However, as discussed in Section 1 and Section 3, it 
should he note that such type of post-processing methods are not variationally consistent and refutes 
the purpose of developing physics-compatible DMP-based meshes. 
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4.4. Sufficient conditions for a rectangular element. For the sake of illustration, consider 
a rectangular element whose vertices are located at (0,0), (o, 0), (a, 6), and (0,6). Our objective is 
to derive conditions on a and 6, such that the local stiffness matrix is weakly diagonally dominant. 
Herein, we consider a pure anisotropic diffusion equation and derive restrictions on the coordinates 
of the rectangular element. Based on the lines of the local stiffness restriction method outlined in 
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Condition No-7. From AM-GM inequality, we have the following relation: 
h Id) aD^i^j 2 . —— 2 1 
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3a ' 36 

which implies that the positive diagonal entries: {Kfj- >0 Vi = 1,2, 3,4, is trivially satisfied. 


Condition No-8. Non-positive off-diagonal entries: {Kfj- <0 d/i j where i = 1,2, 3,4, 
and j = 1,2,3,4, needs to be satisfied. For instance, when i = 1 and j = 2,3,4, this restriction 
gives the following relations: 



< 0 


(4.57a) 


(4.57b) 


6a 2 66 

Additionally, we should also satisfy the restrictions imposed by equations (4.30a)-(4.30c) on dif- 
fusivity tensor. In a similar fashion, one can derive restrictions for other combinations of i and 

j- 


Condition No-9. Weak diagonal dominance of rows: | {Kfj- \ > | {Ke)ij \ '^i,j where 

i = 1,2, 3,4, and j = 1,2,3,4, is trivially satisfied if Condition-7 and Condition-8 are met. This is 
because from equation (4.55), it is evident that i^e)ij = 0 Vi,j where i = 1,2,3,4, 

and j = 1, 2,3,4. 

Finally, as explained in subsubsection 4.2.1, extending the local stiffness restriction approach 
to incorporate advection and linear reaction is straightforward and shall not be dealt with to save 
space. Moreover, it is easy to construct mesh restrictions for any set of arbitrary coordinates of a 
quadrilateral element using symbolic packages like Mathematica. But the resulting inequalities will 
be more complex to analyze mathematically and visualize graphically. 


5. CONCLUDING REMARKS AND OPEN QUESTIONS 

We outlined a general procedure to obtain the restrictions that are needed for a computational 
grid to satisfy various mathematical principles - comparison principles, maximum principles, and 
the non-negative constraint. We illustrated the workings of this procedure by obtaining the mesh 
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restrictions for T3 and Q4 finite elements. The procedure is, however, equally applicable to other 
low-order finite elements. 

First, we critiqued three different approaches to satisfy maximum principles, comparison prin¬ 
ciples, and non-negative constraint for a general linear second-order elliptic equation. A pictorial 
description of a generic relationship between DMPs, DCPs, and NC based on a Venn diagram is 
proposed. This sketch helps to easily relate the space of solutions obtained using mesh restrictions, 
non-negative numerical formulations, and post-processing methods. We then presented necessary 
and sufficient conditions on the stiffness matrix K to meet the mathematical properties. Using 
these conditions, we derived stronger and weaker mesh restrictions for T3 element. The stronger 
mesh restriction corresponds to the anisotropic non-obtuse angle condition while the weaker 
one corresponds to the generalized Delaunay-type angle condition. Motivated by these mesh 
restriction conditions, different kinds of Peclet and Damkohler numbers are proposed for advective- 
diffusive-reactive systems when diffusivity is anisotropic. 

For isotropic diffusivity, we established that acute-angled or right-angled triangle is sufficient to 
satisfy discrete principles. However, for anisotropic diffusivity, we showed that in order to satisfy 
DMPs, DCPs, and NC, all the dihedral angles of a simplex measured in the metric of have 
to be either O (/i||n||cxD,Th + h‘^\\(^\\oo,Th) acute/non-obtuse or O {h\\v\\oo,Th + ^^ll®l|oo,Th) Delaunay. 
Pictorially, this means that the feasible region for T3 and Q4 elements to satisfy various discrete 
principles is based on a metric tensor whose components are a function of anisotropic diffusivity 
with respect to a suitable coordinate system. Then, an anisotropic metric tensor and an iterative 
algorithm to generate various types of DMP-based triangulations are described. Different numerical 
examples and respective DMP-based triangular meshes are presented for different types of D(x) to 
demonstrate the pros and cons of imposing mesh restrictions. Furthermore, the errors incurred in 
satisfying local and global species balance are documented. Based on these numerical experiments, 
the following inferences can be drawn: 

(Cl) For pure isotropic or anisotropic diffusion equation, a coarse DMP-based triangulation is suf¬ 
ficient to satisfy various discrete principles. However, for advection-dominated and reaction- 
dominated scenarios, we need a highly refined DMP-preserving computational mesh to obtain 
non-negative solutions. 

(C2) Existing traditional and non-traditional methods of /i-refinement may not guarantee the sat¬ 
isfaction of DMPs, DCPs, and NC always. 

(C3) On coarse DMP-based meshes, errors incurred in satisfying local and global species balance 
for highly anisotropic diffusion-type problems is considerable due to the skewed nature of 
the mesh elements. Moreover, the decrease in local and global species balance errors upon 
/i-refinement is slow. 

(C4) DMP-based meshes change as one alters the underlying anisotropic diffusivity tensor. 

In the light of the recent developments and motivated by the above discussions, we have chosen 
to emphasize on the following four open problems that we consider particularly interesting in view 
of their mathematical richness, numerical challenges, and potential applications: 

(OPl) In this paper, all the meshes used in the numerical examples are of Delaunay-type. This is 
because most of the existing open source mesh generators such as BAMG, Gmsh, Triangle, 
BL2D, MmgSd, and CGALmesh are Delaunay. Recently, Erten and Ungor [15-17] have 
developed a non-obtuse/acute angled mesh generator called aCute by modifying Triangle. 
However, aCute is restricted to 2D and Eucledian metric tensors. Hence, such a software 
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can only be nsed to satisfy discrete principles for problems involving heterogeneous isotropic 
diffusivity. Having an anisotropic non-obtuse/acute A^-uniform meshing software would be 
of great importance, as the numerical solutions obtained from these meshes not only satisfy 
DMPs, DCPs, and NC, but also converge uniformly (an attractive aspect in finite element 
analysis [12]). To date, there is no such mesh generator. Developing such a software will 
have a profound impact on obtaining physically meaningful numerical solutions for diffusion- 
type equations. 

(OP2) For advection-dominated and reaction-dominated advection-diffusion-reaction problems, mesh 
refinement that adheres to DMPs is needed to obtain stable and sufficiently accurate 
numerical solutions. As described in (C2), not every method of /i-refinement is DMP- 
preserving. Hence, a consistent way of generating a DMP-based /i-refined mesh (that satis¬ 
fies Generalized Delaunay-type angle condition) is still unresolved. 

(OP3) From (C3), it is apparent that local and global mass conservation property is needed. An 
approach to preserve such a property without violating DMPs, DCPs, and NC, is to ob¬ 
tain mesh restrictions for mixed Galerkin formulation based on lowest-order Raviart-Thomas 
spaces. Recently, Huang and Wang [117] have developed a methodology to satisfy DMPs for 
a class of locally conservative weak Galerkin methods using lowest-order Raviart-Thomas 
spaces. However, this methodology is limited to pure anisotropic steady-state diffusion 
equation in two-dimensions. Hence, a mesh restriction based method to satisfy differ¬ 
ent discrete principles, local species balance, and global species balance for anisotropic 
advection-diffusion-reaction equations thus far is unsolved. 

(OP4) In order to construct DMP-based meshes for low-order non-simplicial finite elements such as 
Q4, from subsection 4.4, it is evident that stronger and weaker mesh conditions are needed. 
So far, there are no mesh restriction theorems analogous to simplicial meshes that can 
provide a general framework to construct non-simplicial meshes for anisotropic diffusivity 
tensors. Theoretically and numerically, it would be very interesting and informative to have 
a comparative study on the performance of simplicial vs. non-simplicial DMP-based meshes 
for various benchmark problems discussed in Section 4. 

Nevertheless, due to enormous research activity in the field of advection-diffusion-reaction equa¬ 
tions, it is impossible to list every open question on preserving DCPs, DMPs, and NC. To con¬ 
clude, the research findings in this paper will be invaluable to the research community and finite 
element practitioners in two respects. First, it will guide the existing users on the restrictions to 
be placed on the computational mesh to meet important mathematical properties like maximum 
principles, comparison principles, and the non-negative constraint. Second, for complex geometries 
and highly anisotropic media, this study has clearly shown that placing restrictions on computa¬ 
tional grids may not always be a viable approach to achieve physically meaningful non-negative 
solutions. We hope that this research work will motivate researchers to develop new methodolo¬ 
gies for advective-diffusive-reactive systems that satisfy local and global species balance, comparison 
principles, maximum principles, and the non-negative constraint on coarse general computational 
grids. 
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Figure 2. ABAQUS numerical simulation for an L-shaped domain with multiple 
HOLES: The contours of concentration obtained using ABAQUS are based on three-node 
triangular mesh. 




(a) Non-negative constraint violation (b) Maximum constraint violation 

Figure 3. Minimum and maximum values for concentration in an L-shaped do¬ 
main WITH multiple HOLES: The left and right figures show the minimum and maximum 
values attained in the computational domain based on the numerical simulations for various 
three-node triangular and four-node quadrilateral meshes using ABAQUS. From the above 
hgures, it is evident that the violation in the non-negative and maximum constraints do not 
reduce on mesh refinement. 
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(a) % violation: Non-negative constraint (b) % violation; Maximum constraint 

Figure 4. Percentage of violation in minimum and maximum constraints for 

CONCENTRATION IN AN L-SHAPED DOMAIN WITH MULTIPLE HOLES: The left and right 
hgures show the percentage of nodes that have violated the non-negative and maximum 
constraint for concentration obtained using ABAQUS. These numerical results are based on 
various three-node triangular and four-node quadrilateral meshes. In both of these cases, 
the percentage of mesh nodes that violated these constraints never decreases to zero. 



Figure 5. Venn diagram for the space of solutions based on mesh restrictions: 
A pictorial description of the space of numerical solutions satisfying various DMPs and 
DCPs based on equation (2.25) and Theorem 2.6. It is evident from the above figure that 
VdsmPk C VdsMPr: C VdwmPk C VdwMPk and VdsCPk C VdwCPk- But we would 
like to emphasize that we do not have the following enclosures: VdwCPk C VdwmPk and 
VdsCPk C VdSMPk- 
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Figure 6. Venn diagram for the space of solutions based on various numerical 
FORMULATIONS: A pictorial description of the space of numerical solutions satisfying various 
DMPs, DCPs, and NC. In numerical literature, most of the numerical methods that exist to 
satisfy various DMPs are mainly non-linear. In the past decade, considerable advancements 
have been made to fulfill various version of DMPs for a certain class of linear elliptic and 
parabolic partial differential equations. But it should be noted that there is seldom research 
progress related to satisfaction of different DCPs (see [11, Section 4]). Hence, we would 
like to highlight that developing a general and variationally consistent numerical technique 
to encompass all these discrete principles is still an open problem. 
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Figure 7. Geometrical properties of an arbitrary simplex in 2D: The left figure 
shows a pictorial description of various geometrical properties, such as unit inward normals 
(rip, rig, and rir), dihedral angles in Euclidean metric (/3pg, j3pr, and /3rg), and heights (/ip, 
hq, and hr) of an arbitrary element ftf. £ Th- Correspondingly, the vertices of this triangle 
PQR are given by ip, ig, and Xr- The right figure shows an arbitrary patch of elements 
Sde and Og, (which belong to the triangulation 7h) sharing a common edge Cpg. The edge 
Cpg connects the coordinates ip (= ip') and ig (= ig'). The dihedral angles in Euclidean 
metric opposite to edge Cpg are denoted by /3pg,Oe and (3^^ qi . 




EiGURE 8. T3 ELEMENT FOR HETEROGENEOUS ISOTROPIC DIFFUSIVITY: A pictorial de¬ 
scription of the feasible region (left figure) is shown in light blue color. The right figure 
indicates that the point (a, b) can lie either on the circle with center (5,0) and radius i or 
outside the circular region. The points within the circular region are infeasible. This results 
in two possibilities for choosing a T3 element in the realm of the feasible region, which is 
either a right-angled triangle or an acute-angled triangle. 
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Right-angled triangle 


Acute-angled triangle 


(a) Feasible region for (a, b) 


(b) Possible T3 elements when D^y = 0 


Figure 9. T3 element for anisotropic diffusivity when D^y = 0: A pictorial 
description of the feasible region (left figure) for the coordinates (a, h) is indicated in light 
blue color. The numerical values for the two parameters, which decide the feasible region, 
are chosen to be e = 10 and 77 = 0. In this case, the right Hgure indicates that acute-angled 
and right-angled triangles are possible. As e increases, the coordinate b has to increase 
proportionally to satisfy the inequality given by the equation (4.33a). 


Isotropic 



(d) 


Ri 



Delaunay / non-obtuse/acute 
triangle in Euclidean metric 

R2 



Delaunay-type/non-obtuse/acute 
—1 

triangle in Dq metric 


Figure 10. DMP-based T3 elements for heterogeneous isotropic and 
ANISOTROPIC DIFFUSIVITY: A pictorial description of a mesh generation procedure to ob¬ 
tain a new triangulation using a given background mesh. This new simplicial mesh satishes 
various discrete properties as contrary to the background mesh. The procedure to obtain 
such a triangulation is iterative and is based on Theorems 4.3 and 4.4. 
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(a) Feasible region for (a, b) 



Acute-angled triangle 
Right-angled triangle 
Obtuse-angled triangle 


a 


(b) Possible T3 elements when D^y < 0 


Figure 11. T3 element for anisotropic diffusivity when D^y < 0: The left figure 
indicates the feasible region for the coordinates (a, b) in light blue color. The right figure 
indicates that when D^y < 0, the (Euclidean metric) dihedral angles in the T3 element can 
be acute-angled or right-angled or even obtuse-angled. In this case, we have chosen e = 10 
and r] = —1. For a fixed rj as e increases, the value of coordinate b also increases. So it is 
a daunting task to find a viable T3 element. One can also notice that the feasible region is 
not contiguous. 




Obtuse angled triangle 

Right angled triangle 
Acute angled triangle 


(a) Feasible region for (a, b) 


(b) Possible T3 elements when D^y > 0 


Figure 12. T3 element for anisotropic diffusivity when D^y > 0: The left figure 
indicates the feasible region for the coordinates (a, b) in light blue color. The right figure 
indicates that when Dxy > 0, the (Euclidean metric) dihedral angles in the T3 element can 
be acute-angled or right-angled or even obtuse-angled. In this case, we have chosen e = 10 
and ?7 = 1. For a fixed rj as e increases, the value of coordinate b also increases. For higher 
values of e, it is very difficult to find a suitable T3 element, which can mesh any given 
computational domain. 
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(a) e = 2 and 77 = —1 
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(b) e = 10 and 77 = —1 
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(c) e = 50 and 77 = —1 


(d) e = 100 and 77 = —1 




a 


a 


(e) e = 200 and 77 = —1 


(f) e = 500 and 77 = —1 


Figure 13. T3 element for fixed 77 and varying e: A pictorial description of the 
feasible region (light blue color) for a fixed 77 and varying e. Analysis is performed for 
77 = —1 and e = {2,10,50,100,200,500}. It is evident there is a drastic variation in the 
feasible region as e increases. 
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(a) e = 100 and rj = —2 
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(e) e = 100 and rj = —8 



a 


(b) e = 100 and rj = 2 



a 


(d) e = 100 and rj — 4 



a 


(f) e = 100 and rj = 8 


Figure 14. T3 element for fixed e and varying rj : A pictorial description of the 
feasible region (light blue color) for a fixed e and varying rj . Analysis is performed for 
e = 100 and rj = {—8 ,—4,—2, 2,4,8}. It is evident there is considerable variation in the 
feasible region as rj changes. Also, there is no fixed pattern on this variation about rj . 
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(c) V = ( 0 . 1 , 1 . 0 ) and a = 0 (d) v = ( 0 . 1 , 1 . 0 ) and ol — 1.0 

Figure 15. Test problem The top left figure shows a coarse triangulation (generated 
using Gmsh based on Algorithm 1) employed in the numerical study, which is to the scale. 
The top right figure and the bottom two figures show the concentration profiles obtained 
for various values of the velocity field and linear reaction coefficient using this mesh. 
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60 70 79 89 98 108 118 127 137 146 156 



0.00 0.30 0.60 0.91 1.21 1.51 1.81 2.11 2.42 2.72 3.02 



(a) Element maximum angles: Nv = 539 and Nele = 906 


- 0 . 44 - 0.15 0.14 0.42 0.71 1.00 1.29 1.58 1.86 2.15 2.44 



(b) Delaunay-type condition: v = (0, 0) and a = 0 


- 0 . 43 - 0.15 0.14 0.42 0.71 0.99 1.27 1.56 1.84 2.13 2.41 



(c) Delaunay-type condition: v = (0.1, 1.0) and a = 0 (d) Delaunay-type condition: v = (0.1,1.0) and a — 1.0 

Figure 16. Test problem #1: The top left figure shows the maximum angle possible 
in each element of the mesh. The top right figure and the bottom two figures show the 
element maximum generalized Delaunay-type condition, which is a weaker condition 
as compared to the element maximum anisotropic non-obtuse angle condition. 
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cP(x) = 0.0 




o 
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Figure 17. Test problem #2: The computational domain under consideration is a bi¬ 
unit square with one of its vertices at origin O = (0,0). Homogeneous Dirichlet boundary 
conditions are prescribed on all sides of the square. The volumetric source /(x) is zero 
inside the domain, except for the square region (including the boundaries) located at vertex 
H = (0.375,0.375). In this region, /(x) is equal to unity. Herein, we assume that the velocity 
vector field and linear reaction coefficient are equal to zero everywhere in the computational 
domain. 
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(a) Background mesh: Nv = 47 and Nele = 68 (b) Anisotropic mesh: Nv = 593 and Nele = 1088 

Figure 18. Test problem 4^2-. The left figure shows the background mesh on which 
BAMG operates to give an anisotropic triangulation, which is shown in the right figure. As 
the ratio of the minimum eigenvalue of anisotropic diffusivity tensor to its maximum is 0.1, 
which is not very high, so the resulting triangulation consists of a mixture of skinny and 
normal triangles. 
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Figure 19. Test problem 4^2: The left figure shows the concentration profile based 
on the background mesh, while the right figure shows the concentration profile using the 
anisotropic triangulation. 




(a) Background mesh: Nv = 5079 and Nele = 9918 (b) Anisotropic mesh: Nv = 297 and Nele = 436 

Figure 20. Test problem ^3: The left figure shows the background mesh and the right 
figure shows the anisotropic triangulation obtained using BAMG for all the four cases. 
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(c) Anisotropic mesh: v = (1.5,1.0) and a = 1.0 (d) Anisotropic mesh: v = (5.0, 0.5) and a = 1.0 
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(e) Anisotropic mesh: v = (0, 0) and a — 1000 

Figure 21. Test problem #3: This figure shows the concentration profiles for four 
different cases based on the background mesh and anisotropic meshes shown in Figure 20. 
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(a) Delaunay mesh: Nv = 1564 and Nele = 2826 (b) Delaunay mesh: Nv = 5090 and Nele = 9620 



(c) Delaunay mesh: Nv = 18372 and Nele = 35665 (d) Delaunay mesh: Nv = 69995 and Nele = 137881 

Figure 22. Issues with traditional mesh refinement: Concentration profiles for the 
fracture domain when v = (0.1,1.0) and a = 1.0. The white region in the figures shows the 
area in which the numerical simulation has violated the NC and maximum constraint. 
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(a) Anisotropic mesh: Nv = 8897 and Nele = 17408 (b) Pure anisotropic diffusion: Concentration profile 


Figure 23. Issues with traditional mesh refinement: The left figure shows 
the anisotropic mesh obtained using the traditional mesh refinement procedure on the 
anisotropic triangulation given in Figure 18. The right hgure shows the concentration prohle 
obtained using this refined mesh. It should be noted that this /i-refined DMP-based mesh 
is interiorly connected. Hence, it satisfies DSMPjf. 
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(a) Anisotropic mesh: Nv = 2199 and Nele = 3924 (b) Pure anisotropic diffusion: Concentration profile 

Figure 24. Issues with traditional mesh refinement: The left figure shows 
the anisotropic mesh obtained using the traditional mesh refinement procedure on the 
anisotropic triangulation given in Figure 20. The right hgure shows the concentration prohle 
obtained using this rehned mesh. 
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(a) Anisotropic mesh: Nv = 2647 and Nele = 4789 (b) Pure anisotropic diffusion: Concentration profile 

Figure 25. Issues with non-traditional mesh refinement: The left figure shows a 
refined anisotropic mesh obtained using the non-traditional approach.The right figure shows 
the concentration profile obtained using this refined mesh. It should be pointed out that 
the mesh obtained using this procedure did not converge in Maxiters = 100. Hence, as a 
result, it violates NC and DMPs. 
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(a) Test problem v = (0, 0) and a = 0 
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(b) Test problem ^2\ v = (0,0) and a = 0 
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(c) Test problem ^5^3: v = (0, 0) and a = 0 (d) Test problem #3: v = (1.5, 1.0) and a = 1.0 

Figure 26. Local species balance errors: The figures show the errors incurred in 
satisfying local species balance for various test problems on coarse meshes. 
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